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Abstract 


>*This  report  describes  the  formulation  of  a numerical  procedure  for 
simulating  the  turbulent,  incompressible  flow  over  and  behind  an  axi- 
symmetric,  self-propelled  body.  The  flow  is  treated  in  three  parts: 

(i)  over  the  body,  (ii)  through  the  propeller,  and  (iii)  in  the  very 
near  wake  of  the  body,  where  axial  gradients  cannot  be  ignored. 

Over  most  of  the  body,  the  flow  is  easily  available  from  the  Reynold's 
number  of  the  flow  and  the  body  shape,  but  near  the  body's  tail,  a 
special  treatment  is  required.  The  non-radial  components  of  the  mean 
flow  through  the  propeller  are  obtained  from  a computer  code  based  on  a 
classical  blade  element  analysis,  while  the  radial  component  is  obtained 
analytically  in  terms  of  the  non-radial  components.  The  effect  of  the 
propeller  on  the  aft-body  boundary  layer  turbulence  is  obtained  by 
integrating  the  very  neaT  wake  turbulence  model  equations  across  the 
plane  of  the  propeller. 

In  the  very  near  wake,  the  flow  is  computed  from  the  time-averaged 
Navier-Stokes  equations,  closed  with  a second-order  model  for  the  turbulence 
correlations.  This  system  of  equations  is  elliptic,  so  conditions  at 
all  boundaries  of  the  very  near  wake  region  are  required.  In  order  to 
choose  a sufficiently  large  solution  domain  without  sacrificing  adequate 
resolution  just  behind  the  propeller,  the  equations  are  transformed 
logarithmically.,.  A uniformly  discretized  version  of  the  transformed 
equations  is  solved  using  an  ADI  (Alternating-Direction- Implicit)  procedure. 
The  flow  is  treated  \as  transient  and  is  allowed  to  evolve  to  a steady 
state.  \ 

Three  sample  calculations  are  presented  with  comparisons  to  experi- 
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coordinates  and  symbols  appropriate  for  boundary  layer  flows.  Section  3 
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1 . Introduction 

When  modeling  the  flow  around  and  behind  a self-propelled 
body,  it  is  convenient  to  divide  the  domain  of  the  flow  Into  two  regions 
(see  Figure  1.1).  These  regions  are  Region  I,  which  contains  the  flow 
around  the  body  and  in  its  very  near  wake,  where  axial  gradients  of  flow 
quantities  are  important,  and  Region  II,  the  downstream  portion  of  the 
wake,  where  axial  gradients  can  be  neglected.  This  division  has  important 
economic  consequences,  since  an  expensive  detailed  calculation  of  the 
flow  is  needed  only  in  Region  I,  the  results  of  which  can  be  used  to 
provide  a plane  of  initial  conditions  for  a relatively  Inexpensive, 
marching  type  calculation  in  Region  II. 

The  phenomena  associated  with  the  boundary- layer,  propeller, 

and  near-wake  portions  of  Region  I (what  we  have  termed  region  14  in 

Figure  1.1)  have,  until  now,  received  little  attention,  but  accurate 

flow  calculations  in  this  region  are  essential  in  order  to  correctly 

initialize  calculations  in  Region  II.  It  is  the  purpose  of  this  report 

to  describe  a method  for  modeling  the  flow  in  Region  1*5,  the  key  elements 

of  which  include  (Figure  1.2):  (1)  the  boundary-layer  flow  near  the 

stern  of  the  body  (z  < 0),  (2)  the  flow  through  the  propeller 

(0  < z < 0+),  and  (3)  the  flow  in  the  near-wake  region  (0+  < z < zq). 

Here,  z is  the  distance  downstream  of  the  body's  tail  and  z marks  the 

* ° 
beginning  of  Region  II.  The  above  flow  elements  are  treated  in  Sections 

2.1,  2.2,  and  Section  3,  respectively,  of  this  report.  Before  proceeding 

to  these  discussions,  we  first  briefly  summarize  the  characteristics  of 

each  element. 

The  body  boundary-layer  flow  is,  in  general,  easily  computed 
or  estimated  from  the  flow  Reynolds  number  and  body  shape.  Near  the 
stern  of  the  body,  however,  where  the  body  radius  becomes  small  compared 
with  the  turbulent  boundary- layer  thickness,  the  flow  evaluation  is  not 
straight-forward  and  a special  treatment  is  required.  As  the  flow 
passes  through  the  propeller,  both  the  mean  and  the  turbulent  flow 
fields  are  altered.  The  flow  at  z - 0+  , required  for  the  near-wake 
calculation,  is  evaluated  in  three  parts:  (a)  the  mean  axial  and  swirl 

*A  typical  value  for  zq  , based  on  Flow  Research  Laboratory  results 
(Lin,  et.  al  (1974)),  is  12  body  radii  behind  the  body  tail. 


* " 
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velocities,  (b)  the  mean  radial  velocity,  and  (c)  the  Reynolds  stresses. 
In  the  near-wake  region,  the  flow  is  characterized  by  relatively  large 
radial  and  axial  mean  shears  and  high  energy  exchange  rates;  it  is 
strongly  influenced  by  the  body  and  propeller  configurations.  In 
particular,  significant  production  of  turbulent  energy  arises  from  the 
swirl  induced  velocity  gradients  and  the  non-uniformities  in  the  axial 
flow.  Accompanying  the  increasing  turbulent  energy  is  an  increase  in 
the  rate  at  which  turbulent  energy  is  dissipated  into  heat.  In  fact, 
according  to  the  measurements  of  Gran  (1972),  the  bulk  of  the  excess 
energy  input  to  the  flow  field  by  the  propeller  is  dissipated  in  the 
near-wake  region. 

A computer  code,  ICWAKE,  based  on  the  methods  described  in 
this  report  has  been  written  to  compute  the  flow  in  Region  lh  and  thus 
provide  initial  conditions  for  the  Region  II  calculation.  The  code  is 
described  by  Grabowskl  and  Robins  (1976) . 
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2.  Propeller  Region 

2.1  Aft-Body  Boundary-Layer  Characterization 

In  order  to  estimate  the  fluid  flow  velocities  in  the  near 
wake,  consideration  must  first  be  given  to  the  boundary- layer  flow  field 
which  envelopes  the  surface  of  the  body.  The  body  shape  is  assumed  to 
be  axisymmetric  and  slender  in  the  sense  that  the  body  length  is  considerably 
larger  than  its  maximum  frontal  diameter.  A typical  characterization  of 
the  body  shape  is  given  by  a Rankine  Ovoid  with  a length-to-diameter 
ratio  of  about  10.  Often,  the  aft-end  of  the  body  is  more  gradually 
tapered  than  a Rankine  Ovoid  in  order  to  alleviate  the  tendency  for  the 
boundary  layer  to  separate.  The  orientation  of  the  major  axis  of  the 
body  is  assumed  to  coincide  with  the  direction  of  the  approaching  flow. 

Slight  deviations  (i.e.,  small  angles  of  attack)  are  not  expected  to 
grossly  affect  the  wake  flow  field.  , 

The  boundary  layer  which  develops  on  the  body  is  assumed  to  be 
fully  turbulent  well  ahead  of  the  mid-section  which  implies  that  the 
Reynolds  number  (based  on  the  speed  of  the  body  and  its  frontal  diameter) 
is  in  excess  of  10^.  Over  most  of  the  body  surface,  the  description  of 
the  boundary  layer  is  amenable  to  conventional  analysis  because  it  is 
usually  thin  in  comparison  with  the  local  body  radius  (6/R«l).  Near 
the  stern  of  the  body,  its  description  becomes  difficult.  This  difficulty 
arises  because  of  two  separate  but  concurrent  effects:  because  of  the 
decreasing  cross-section  body  radius,  the  boundary  layer  undergoes  a 
change  from  a flat-plate-like  flow  (6/R«l)  to  one  where  axisymmetric 
effects  dominate  and  6/R-*»  ; and  because  there  exists  substantial 
interaction  between  the  thick  boundary  layer  near  the  aft  end  of  the 
body  and  the  potential  flow  outside  it,  the  boundary  layer  is  subjected 
to  both  streamwise  and  transverse  pressure  gradients.  Thus,  several 
assumptions  required  by  conventional  boundary-layer  theory  are  violated. 
Despite  the  considerable  amount  of  effort  which  has  been  devoted  to  the 
development  of  an  analytical  model  to  describe  the  aft-body  boundary- 
layer  characteristics  (e.g.,  Patel  1973a  and  1973b),  these  tools  have 
not  yet  been  shown  to  be  reliable  for  practical  use  at  very  high  Reynolds 
numbers . 


» 
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In  lieu  of  employing  unproven  computer  codes  to  calculate  the 
aft-body  boundary  layer  flow,  we  have  examined  the  available  experimental 
data  and  employed  somewhat  more  conventional  extrapolation  procedures  in 
order  to  estimate  the  boundary-layer  characteristics  of  a body  Immersed 
in  a very  high  Reynolds  number  flow.  In  the  following  we  first  examine 
the  mean  flow  field  characteristics  and,  then,  examine  the  turbulence 
characteristics . 

2.1.1  Mean  Flow  Field 

In  our  first  attempt  at  correlating  the  mean  velocity  profile 

of  the  aft-body  boundary  layer. we  employed  the  conventional  boundary- 

w w ^UT 

layer_correlation  schemes  such  as  plotting  — (y/6)  , — (— — ) 

and  -^e  (y/6)  where  w is  the  streamwise  fluid  speed,  is  the 

local  freestream  speed,  uT  is  the  friction  velocity,  y is  the  distance 

measured  normal  to  and  from  the  body  surface  and  6 is  the  local 

boundary- layer  thickness.  Such  graphs  (using  Patel's  1974  experimental 

data)  failed  to  display  any  noticeable  degree  of  collapse  to  a single 

curve.  The  lack  of  collapse  can  be  attributed  to  the  boundary-layer/ 

mean-flow  interaction  and  to  the  failure  of  the  flow  to  satisfy  the 

boundary-layer  assumption  6/R«l  at  the  trailing  edge. 

We  then  hypothesized  that,  in  spite  of  the  complex  dynamics  of 

the  flow  in  the  aft-body  region,  the  mechanisms  whereby  the  total  mean 

flow  energy  is  converted  to  turbulent  energy  might  not  be  drastically 

different  from  the  mechanisms  which  act  in  a flat-plate  boundary  layer. 

To  test  this  hypothesis,  Patel's  (1974)  very  extensive  measurements  of 

the  aft-body  boundary  layer  were  used  to  compute  the  total  mean  flow 

2 

energy  (i.e.,  (p/p)+(q  /2)  , where  the  pressure  term,  p/p  , represents 

2 

the  potential  energy  of  the  fluid  and  the  velocity  term  q /2  represents 
the  kinetic  energy  of  the  fluid)  for  each  streamline  at  a number  of 
measurement  stations  along  the  surface  of  the  body  which  he  utilized. 

A graph  of  this  reduction  of  Patel's  data  is  shown  in  Figure  2.1  where 
the  ordinate  is  the  square  root  of  the  stagnation  pressure  coefficient 

p-poo  /q  \2 

(C  ■ =■  + ( — I ) , and  the  abscissa  Is  the  streamfunctlon  , 

p Jipw l \W»/ 
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normalized  by  streamf unction,  i|>e  , at  the  "edge"  of  the  boundary  layer 
(where  Cp  - 0.99).  Here  the  subscript  <*>  denotes  free-stream  values. 
The  value  of  ^ at  any  survey  location  was  determined  from  Patel's  data 
from 


y 

< l>  m f 2ir(y'+y  )q*n  dy' 

0 D 

where  y^  is  the  distance  from  the  body  axis  to  the  body  surface  measured 
along  the  local  normal  and  n is  a unit  vector  parallel  to  the  body 
surface  at  the  point  where  the  normal  intersects  it. 

The  rationale  for  this  particular  choice  of  coordinate  scales 
will  be  discussed  below.  It  is  evident  from  Figure  2.1  that  all  of  the 
computed  values  lie  within  a very  narrow  band  in  spite  of  the  fact  that 
the  local  pressure  gradient  and  ratio  of  boundary-layer  thickness  to 
body  radius  both  change  markedly  over  the  range  of  survey  stations  (the 
numbers  shown  next  to  the  tabulated  symbols  represent  the  fraction  of 
the  body  length  measured  from  the  leading  edge  where  the  boundary- layer 
survey  was  conducted).  At  z*/L  * 0.662  , for  example,  the  streamwise 
pressure  gradient  is  negligible  and  6/R  *0.15  , whereas  beyond  z*/L  * 0. 
the  boundary  layer  is  subject  to  an  adverse  pressure  gradient  and  at 
z*/L  « 0.99  , 6/R  * 13  . Here,  L is  the  body  length  and  z*  is  the 

distance  behind  the  body  nose. 

In  order  to  lend  additional  credence  to  this  correlation,  the 
universal,  flat-plate,  boundary- layer  profile  given  by  Coles  and  Hirst 
(1969)  is  shown  by  the  heavy  line  in  Figure  2.1  after  converting  it  to 
the  coordinates  used  in  the  figure.*  The  close  correspondence  between 
Coles'  flat-plate  profile  and  the  aft-body  boundary-layer  measurements 
taken  by  Patel  supports  the  original  hypothesis.  Two  examples  of  one 
other  flat-plate  boundary  layer  correlation  scheme  are  also  shown  in 
Figure  2.1;  these  examples  are  for  the  case  where  the  velocity  profile 
is  approximated  by 


A representative  value  of  the  friction  coefficient  Cf  measured  by  Patel 
has  been  used  to  compute  this  curve;  it  should  be  noted  that  the  measured 
Cf  is  some  30%  lower  than  would  be  estimated  for  a flat-plate  boundary 
layer  at  a Reynolds  number  based  on  the  body  half-length.  This  deviation 
is  most  likely  a result  of  the  boundary  layer  being  tripped  near  the 
leading  edge. 
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w/wg  * (y/6)1^11  . 

Profiles  of  this  representation  for  n = 5 and  7 are  indicated 
by  the  dashed  lines  when  computed  in  U.e  coordinates  shown  in  Figure  2.1. 
Although  this  representation  has  less  theoretical  basis  than  Coles',  the 
case  for  n = 5 appears  to  fit  the  data  quite  well,  at  least  for  the 
inner  50%  of  the  boundary-layer  mass  flux.  The  fact  that  the  l/7th 
power  law  profile  appears  to  represent  the  data  somewhat  better  in  the 
outer  region  again  may  be  attributed  to  the  boundary-layer  trips. 

The  comparison  between  these  analytic  velocity  profiles  and 
the  experimental  data  gives  an  indication  as  to  why  the  particular 
coordinates  were  chosen  to  plot  the  data  and  the  curves.  For  a flat- 
plate  boundary-layer,  the  abscissa  represents  nothing  more  than  the 
velocity  ratio  within  the  boundary  layer  since  p = p^  . This  would  be 
true  also  of  the  aft-body  boundary  layer  if  the  flow  were  to  be  inviscid 
and  were  to  turn  parallel  to  the  body  axis  with  the  static  pressure 
recovering  to  the  free  stream  pressure.  If  one  desires  the  boundary- 
layer  velocity  profile  in  the  aft-body  region,  it  is  necessary  that  the 
pressure  distribution  be  known  ji  priori.  Often  this  is  not  known, 
however,  because  of  the  influence  of  some  other  phenomenon  such  as  a 
marine  propeller  mounted  at  the  trailing  edge.  Even  if  this  could  be 
estimated  by  some  technique,  the  correlation  shown  in  Figure  2.1  would 
only  permit  the  calculation  of  the  magnitude  of  the  fluid  velocity 
within  the  boundary  layer;  the  flow  direction  could  not  be  computed. 

It  should  be  added  that  some  of  the  other  data  obtained  by 
Patel,  et  al.,  indicate  that  the  flow  direction  within  the  boundary 
layer  varies  linearly,  from  being  parallel  to  the  surface  at  the  body  to 
being  parallel  to  the  undisturbed  flow  at  the  boundary  layer  edge. 

The  normalization  of  the  abscissa  by  i|>e  simply  accounts  for 
the  fact  that  even  a flat-plate  boundary  layer  entrains  fluid  as  it 
develops  downstream.  Manual  calculations  (not  shown  here)  of  the 
entrainment  rate  for  a flat-plate  boundary  layer  and  for  an  aft-body 
boundary  layer  indicated  that  each  entrains  fluid  at  about  the  same  rate 
in  spite  of  the  geometric  and  pressure  gradient  differences  between 
them.  Thus,  using  the  flat-plate  computation  of  the  amount  of  entrained 


I 
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fluid,  the  quantity  can  be  estimated  at  any  location  along  the  body 

surface  from  (see  Appendix  A.l) 

\pe(z*)  * 4 wg  A(z*)(Cp)F 

where:  A(z*)  is  the  wetted  surface  area  of  the  body  between  the 

leading  edge  and  the  location  z*  . 

(Cp)F  p is  the  flat-plate  value  of  the  total  friction 

coefficient  computed  using  the  Reynolds  number  w z*/6 

e 

(possibly  modified  in  the  event  a boundary layer  trip  is 
employed) . 


r * — ?ft- 


The  line  of  reasoning  for  the  apparent  independence  of  the 
entrainment  rate  on  surface  geometry  is  that  it  is  roughly  proportional 
to  the  turbulent  intensity.  For  a flat-plate  boundary- layer  it  has  also 
been  observed  that  the  turbulent  intensity  is  about  three  times  the 
local  friction  velocity  (u^  = wg  Jc^/2) . For  the  aft-body  flow,  however, 
Patel's  measurements  indicate  that  this  latter  relation  is  no  longer 
valid  but  that  the  turbulent  intensity  remains  at  about  the  same  level 
even  though  Cp  is  considerably  smaller. 

In  spite  of  the  apparent  success  of  the  correlation  of  the 

aft-body  boundary-layer  data  with  flat-plate  predictions,  the  calculation 

of  the  physical  velocity  profile  (i.e.,  w/w  (r))  is  not  a trivial 

e 

matter  even  if  the  turning  of  the  flow  parallel  to  the  axis  and  the 
recovery  of  the  static  pressure  (with  viscous  effects  neglected)  could 
be  approximated  (this  turning  and  pressure  recovery  is  exactly  what  is 
usually  assumed  for  computing  the  flow  through  a marine  propeller  such 
as  that  given  in  Section  2.2).  In  this  event,  the  streamline  location 
can  be  computed  from 
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and,  by  cross-plotting,  the  velocity  profile  w/wfi(r)  can  then  be 
determined.  If  a power-law  boundary- layer  velocity  profile  correlation 
were  adequate  so  that 
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then  at  the  trailing  edge  of  the  body,  the  above  Integral  may  be  evaluated 
with  the  result  that 


Note  that  this  implies  that  the  velocity  profile  when  expressed  as  a 
function  of  r/6  is  less  full  than  a flat-plate  profile  expressed  as  a 
function  of  y/6  . 

2.1.2  Turbulence  Characteristics 

The  data  of  Patel,  et  al.  (1974)  have  also  been  examined  to 
estimate  the  characteristics  of  the  turbulence  within  the  aft-body 
boundary  layer.  This  effort  has  involved  converting  Patel's  measure- 
ments into  the  form  of  the  dependent  variables  which  are  used  in  the 
ICWAKE  code.  The  reader  is  referred  to  Section  3 for  the  precise  defini- 
tion of  these  dependent  variables.  For  the  purposes  here,  it  is  noted 

that  R represents  the  expected  value  of  the  mean-square  velocity 
z z 

fluctuation  in  the  direction  parallel  to  the  body  axis.  The  quantities 

Rq0  , Rrr  , Rrz  , and  the  other3  are  defined  similarily. 

In  our  effort  to  correlate  these  data,  we  have  taken  the 

liberty  to  incorporate  some  Reynolds  number  scaling  into  the  graphs 

discussed  below.  This  involves  nothing  more  than  incorporating  the 

approximate  Reynolds  number  scaling  for  flat-plate  boundary  layers  where 

all  turbulent  stresses  are  proportional  to  the  square  of  the  flat-plate 

friction  velocity  into  the  aft-body  boundary-layer  correlations.  That 

2 

is,  in  a flat-plate  boundary  layer  It  is  found  that  R /u  is  a function 
of  y/6. 

Based  on  the  reasoning  presented  in  Section  2.1.1  regarding 
the  turbulent  entrainment  rate,  we  assume  here  that  even  in  the  aft-body 
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boundary  layer  the  magnitudes  of  the  turbulent  stresses  are  proportional 

to  the  flat-plate  friction  velocity  computed  at  the  Reynolds  number 

w<x>z*/v  . The  justification  for  this  assumption  is  that  it  works  for 

flat-plate  boundary  layers  and  that  it  agrees  fairly  well  with  the 

turbulent  entrainment  rate  for  the  aft-body  boundary  layer.  It  should 

be  noted  that  this  Reynolds  number  scaling  only  introduces  a slight 

2 

correction  to  the  stresses  since  u_  - (C-)_  D and  the  local  flat- 

T t fir* 

plate  friction  coefficient,  C , only  changes  a factor  of  two  or  so 

r 

over  several  orders  of  magnitude  change  in  Reynolds  number. 

The  distribution  of  the  turbulent  intensity  across  the  aft- 
body  boundary  layer  as  measured  by  Patel,  et  al.,  is  shown  in  Figure 
2.2.  On  this  figure  the  turbulent  intensity  (defined  as  the  square  root 
of  one-half  the  sum  of  the  three  turbulent  normal  stresses)  is  plotted 
as  a function  of  the  fractional  distance  across  the  local  boundary-layer 
thickness  for  the  four  survey  stations  reported  by  Patel  (and  represented 
by  the  open  symbols).  Klebancff's  (1955)  flat-plate  measurements  which 
are  represented  by  the  solid  circles  are  also  shown  in  Figure  2.2. 

Except  for  the  region  near  the  body  surface  ((y-yj))/6  ~ 0.2)  at  the 
trailing  edge,  these  data  all  lie  within  a band  whose  mean  is  indicated 
by  the  solid  line  shown  in  Figure  2.2. 

The  degree  of  anisotropy  of  the  aft-body  boundary-layer 
turbulence  is  shown  in  Figure  2.3  where  the  ratio  of  the  radial  to  the 
axial  and  the  ratio  of  the  circumferential  to  axial  turbulent  intensities 
are  shown  as  a function  of  position  within  the  boundary  layer.  Although 
there  is  considerable  scatter  in  both  of  these  sets  of  data,  the  crude 
approximation  that  both  the  radial  and  the  circumferential  intensities 
are  about  60  percent  of  the  axial  turbulent  Intensity  over  most  of  the 
aft-body  boundary  layer  is  not  too  unreasonable. 

The  distribution  of  turbulent  shear  stress  R across  the 

rz 

aft-body  boundary  layer  is  shown  in  Figure  2.4.  After  several  other 
correlation  attempts  (such  as  estimating  the  eddy  viscosity  rather  than 
the  turbulent  stress),  it  was  our  opinion  that  Bradshaw's  (1967)  normali- 
zation, in  which  the  stress  is  assumed  proportional  to  the  turbulent 
energy,  appeared  to  give  the  best  correlation  of  Patel's  data.  It 
should  be  observed  that  for  the  data  at  z*/L  - 0.662  (where  the 
boundary  layer  is  thin  and  the  pressure  gradient  is  small),  Patel 'a 
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data  agree  quite  well  with  Bradshaw's  estimate.  As  the  boundary  layer 

thickens  and  decelerates  toward  the  aft-end  of  the  body,  the  shear 

stress  in  the  outer  portion  of  the  boundary  layer  decreases  relative  to 

Bradshaw's  flat-plate  boundary-layer  correlation.  Patel's  turbulent 

shear  stress  data  appear  to  be  reasonably  correlated  by  the  heavy  dashed 

line  shown  in  Figure  2.4.  A Reynolds  number  scaling  of  the  shear  stress 

data  shown  in  Figure  2.4  is  implicit  in  that  the  shear  stress  normali- 
2 2 

zation  quantity  k is  scaled  with  uT  (Figure  2.2). 

The  other  two  turbulent  shear  stresses  R « and  R » are 

r6  z9 

zero  within  the  boundary  layer.  It  should  be  remarked  that  if  the 
turbulence  passes  through  a marine  propeller  (which  exerts  a step  change 
in  the  swirl  velocity  at  the  plane  of  the  propeller),  these  cwo  stress 
components  will  no  longer  be  identically  zero.  The  effect  of  a propeller 
on  all  of  the  turbulent  velocity  correlations  is  discussed  further  in 
Section  2.2.2. 

By  far  the  most  difficult  correlation  that  was  to  be  extracted 
from  Patel's  data  was  for  the  turbulent  dissipation  rate.  The  initial 
distribution  of  this  quantity  aft  of  the  propeller  is  required  by  the 
ICWAKE  code  in  Section  3.  This  quantity  was  not  measured  nor  could  it 
be  inferred  directly  from,  say,  the  decay  rate  of  the  turbulent  intensity 
along  a given  streamline.  To  proceed,  we  adopted  a model  for  the  local 
turbulent  dissipation  rate,  such  as  given  by  Tennekes  and  Lumley  (1973). 
Basically,  this  assumes  that  the  dissipation  rate  £ 1b 


e 


/A 


where  A is  the  longitudinal  Integral  scale  size  of  the  turbulent 

velocity  fluctuations.  In  the  aft-body  region  we  have  assumed  that  A 

is  no  longer  proportional  to  the  actual  boundary- layer  thickness  but, 

rather,  it  is  proportional  to  the  thickness  of  an  undlatorted  boundary 

layer.  Furthermore,  the  quantity  R has  been  assumed  to  be  proportional 
2 22 

to  (u  )_  according  to  the  correlation  shown  in  Figure  2.2.  With 

this  line  of  reasoning,  one  is  then  led  to  the  conjecture  that 
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To  test  this  contention,  Patel's  data  were  used  to  compute  the 

3/2 

left-hand-side  of  the  equation  (which  is  nearly  k times  a very  weak 

function  of  Reynolds  number  to  account  for  the  slight  change  in  (u  )_ 

T r »r  • 

and  6 _ for  the  different  survey  stations  examined  by  Patel).  The 

r • * • 

results  of  this  computation  are  shown  in  Figure  2.5.  Note  that  the 
scatter  of  the  data  is  considerable  and  that  it  is  also  somewhat  different 
than  the  data  for  a flat-plate  boundary  layer.  In  spite  of  this,  we 
have  attempted  to  estimate  the  turbulent  dissipation  rate  distribution 
across  the  aft-body  boundary  layer  in  a manner  indicated  by  the  heavier 
dashed  line.  Obviously,  an  improved  correlation  would  require  not  only 
a better  underlying  scheme  but  more  detailed  and  varied  measurements. 

Whether  a better  correlation  than  that  given  above  is  really  needed  or 
not  depends  on  how  sensitive  the  wake  calculations  are  with  respect  to 
this  quantity. 

2 . 2 Influence  of  a Marine  Propeller  on  the  Aft-Body  Boundary  Layer 

In  the  previous  subsection  of  this  report,  the  characteristics 
of  the  aft-body  boundary  layer  were  examined  and  correlated  so  that  one 
could  make  reasonable  estimates  of  the  mean  velocity  field  and  turbulence 
within  the  fluid  prior  to  the  formation  of  the  wake.  This  correlation 
scheme,  it  should  be  noted,  is  extracted  from  data  obtained  for  the  flow 
field  of  a body  whose  speed  is  maintained  by  towing  it  through  the 
surrounding  fluid.  For  most  cases  of  practical  interest,  however,  the 
body  propels  Itself  through  the  surrounding  fluid  by  means  of  an  aft- 
mounted  marine  propeller.  Because  of  the  presence  of  a propeller  one 
might  question  the  applicability  of  the  towed-body,  boundary-layer  data 
to  a more  realistic  case. 

Experimentally,  it  has  been  observed  that  the  upstream  influence 
of  a marine  propeller  is  limited  to  a distance  on  the  order  of  the 
propeller  diameter.  This  influence  is  judged  by  the  change  in  the 
surface  pressure  distribution  along  the  body  operating  with  and  without 
a propeller  (see  McLemore,.  1962,  for  example).  Thus,  the  development  of 
the  body's  boundary  layer  is  unaffected  by  the  presence  of  a propeller 
upstream  of  this  region  which  occupies  only  the  last  5 percent  or  so  of  the 
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vehicle's  length.  The  mean  flow  field  and  turbulence  characteristics 
described  in  the  previous  section  are  presumably  valid  over  most  of  the 
vehicle  surface.  Within  this  latter  5 percent  of  the  body  length,  the 
influence  of  the  propeller  is  such  that  the  fluid  is  accelerated  towards 
the  propeller  because  of  the  reduced  upstream  pressure ■ field  which 
accompanies  the  propeller  thrust.  A more  detailed  description  of  this 
upstream  influence  of  the  propeller  will  be  given  in  Section  2.2.1. 
Qualitatively,  however,  the  description  here  implies  that  only  the  mean 
flow  field  is  affected  by  the  propeller;  the  turbulence  within  the  aft- 
body  boundary  layer  is  assumed  to  be  simply  convected  along  with  the 
mean  flow.  Furthermore,  because  the  upstream  influence  of  the  propeller 
is  so  limited  in  extent,  it  is  not  expected  that  the  energy  exchange 
processes  associated  with  the  boundary- layer  turbulence  will  have  sufficient 
time  to  exert  a significant  influence  on  the  mean  flow.  Hence,  we 
assume  that  the  entire  flow  immediately  upstream  of  the  propeller  is 
inviscid  (but  rotational)  and  that  each  turbulence  correlation  can  be 
treated  as  a passive  scalar.  Immediately  downstream  of  the  propeller  an 
assumption  like  this  would  be  invalid  since  comparatively  large  velocity 
gradients  exist  because  of  the  propeller  generated  axial  and  swirl 
velocities  (the  latter  of  which  does  not  exist  upstream) . These  mean 
flow  velocity  gradients  downstream  of  the  propeller  increase  the  rate  of 
energy  exchange  from  the  mean  flow  to  the  turbulence  thus  increasing  its 
intensity  and  the  influence  of  the  dissipative  and  diffusive  effects  of 
the  turbulence.  This  is  why  a complicated  turbulence  model  computer 
code  is  necessary  in  the  near-wake  but  is  unnecessary  upstream  of  the 
propeller. 

2.2.1  Propeller  Induced  Mean  Flow  Field 

Regarding  the  Influence  of  a marine  propeller  on  the  rotational 
(but  assumed  inviscid)  flow  entering  the  propeller,  a number  of  analytical 
investigations  have  been  made.  In  particular,  the  analysis  developed  by 
Schwartz  and  Bernstein  (1975)  was  examined.  Their  model  follows  the 
classical  propeller  analysis  of  combined  blade  element,  momentum  balance, 
and  vortex  theory.  A computer  code  (op.  clt.)  based  on  their  model 
predicts  the  Induced  axial  and  swirl  mean  velocities  at  the  propeller 
plane  for  a given  propeller  (size,  number  and  shape  of  blades,  rotational 
speed)  and  a given  undisturbed  upstream  boundary-layer  velocity  profile. 
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Two-dimensional  airfoil  theory  (with  corrections  for  finite  span)  is 
used  for  calculating  the  induced  velocities  within  constant  radius 
annuli.  The  mean  radial  velocity  at  the  propeller  plane  is  not  de- 
termined by  the  Schwartz-Bernstein  model. 

An  alternate  derivation  of  the  upstream  influence  of  the 
propeller  is  given  in  Appendix  A. 2.  This  analysis  leads  to  an  exact 
solution  for  the  nearly  inviscid  linearized  equations  describing  the 
flow  through  a marine  propeller.  A significant  outcome  of  this  analysis 
is  that  it  yields  an  expression  for  the  mean  radial  velocity  at  the 
propeller  plane  in  terms  of  the  induced  axial  velocity  there.  This  is 
an  Important  result  because  it  turns  out  (Appendix  A. 2),  within  the 
approximations  of  the  theory , that  the  radial  flow  speed  is  of  the  same 
order  of  magnitude  as  the  induced  axial  flow  speed,  and  hence  cannot  be 
ignored.  (Note  that  the  induced  axial  flow  speed  is  the  change  in  axial 
speed  from  "far"  upstream  of  the  propeller  to  the  propeller  plane, 
which,  depending  on  the  propeller  and  its  thrust,  can  be  as  much  as  20 
percent  of  the  freestream  speed.) 

On  the  other  hand,  it  was  determined  that  even  though  the  new 
analysis  is  in  some  ways  more  complete  than  the  previous  work,  the  new 
estimates  for  the  induced  axial  and  swirl  velocities  do  not  differ 
greatly  from  the  previous  estimates.  Thus,  the  final  computer  model  for 
the  mean  flow  in  the  propeller  plane  was  chosen  to  be  the  Schwartz- 
Bernstein  code  with  the  addition  of  a routine  to  compute  the  mean  radial 
flow  speed.  The  equation  that  relates  the  radial  flow  speed  to  the 
axial  flow  speed  change  is  (Appendix  A. 2) 

00 

u(r)  - - J Aw(p)Gu(r,p)pdp  , 
o 

where  Aw(p)  is  the  Induced  axial  velocity  at  the  plane  of  the 
propeller  as  a function  of  radial  distance,  and  is 
obtained  from  the  Schwartz~Bernstein  model. 


and  Gu(r,p)  la  a "Green's  function"  given  by 


G (r  p)  .L-  [E&L  KimJ_] 
uU,p;  nr  [r  - p r + pj  » 


where  K(m)  and  E(m)  are  complete  elliptic  functions  of  the  1st  and 
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2nd  kind,  respectively,  and  the  modulus  m = 4 rp/(r  + p)  . Examination 
of  the  function  shows  that  an  integrable  singularity  exists  at 

P - r . 

2.2.2  Effect  of  Propeller  on  Aft-Body  Boundary-Layer  Turbulence 

As  explained  in  the  previous  two  sections  of  this  report,  the 
turbulence  characteristics  of  the  aft-body  boundary  layer  are  not 
expected  to  change  substantially  in  the  region  upstream  of  the  propeller. 
That  is,  the  turbulence  is  assumed  to  be  simply  convected  along  with 
the  flow  through  this  region  until  it  arrives  at  the  propeller  plane. 

As  the  turbulence  passes  through  the  propeller,  however,  it  could  be 
markedly  changed  by  the  mechanical  "chopping"  of  the  blades;  also  the 
mean  propeller  flow  field  could  substantially  distort  the  initial  wake 
turbulence  and  some  turbulence  could  be  added  to  the  body  wake  from  the 
boundary  layers  on  the  propeller  blades.  This  is  particularly  important 
if  the  propeller  is  operating  in  a regime  whereby  the  blade  boundary 
layers  separate  from  the  blade  surface.  It  is  the  purpose  of  this 
section  to  examine  the  potential  effect  on  the  initial  wake  turbulence 
level  of  each  of  these  three  mechanisms;  mechanical  chopping,  distortioi 
of  the  incoming  turbulence,  and  generation  of  new  turbulence. 

Mechanical  Chopping 

As  a turbulent  eddy  is  convected  through  a marine  propeller  by 
the  mean  flow,  it  may  be  cut  into  two  or  more  segments  by  a blade.  Because 
the  propeller  blades  are  lifting  surfaces,  the  fluid  which  passes  over 
the  suction  side  of  the  blade  (upstream  side)  does  so  in  less  time  than 
that  which  passes  on  the  pressure  side.  Therefore,  downstream  of  the 
propeller  the  segments  of  the  eddy  do  not  rejoin  but  are  displaced 
relative  to  one  another  by  a distance  note  where  a is  the  blade 
angle-of-attack  and  c is  the  blade  chord  (this  estimate  is  easily  made 
using  conventional  airfoil  theory).  Further,  the  axial  distance  between 

1 W» 

the  passage  of  two  adjacent  blades  is  jj  where  B is  the  number  of 

blades,  and  n is  the  propeller  rotational  speed  (in  cycles  per  unit 

time).  If  typical  values  for  the  number  of  blades,  the  blade  chord 

length,  etc.  are  inserted  into  these  relations,  it  may  be  estimated  that 

the  net  displacement  of  a severed  turbulent  eddy  is  on  the  order  of  0.15 

D and  the  distance  between  blade  passages  is  about  0.1  D where  D 
P P P 
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1s  the  diameter  of  the  propeller.  For  eddy  sizes  much  larger  than  0.1 
Dp  the  splitting  of  the  eddy  is  negligible  whereas  for  eddy  sizes 
much  less  than  0.1  Dp  only  a small  fraction  will  experience  a permanent 
split.  It  is  concluded  that  mechanical  chopping,  at  best,  only  causes  a 
very  localized  redistribution  of  the  turbulent  energy  spectrum  in  the 
neighborhood  of  those  eddy  sizes  which  are  about  0.1  D . This  effect 
can,  therefore,  be  ignored  for  the  present  purposes. 

Generation  of  Turbulence  Within  the  Propeller  Blade  Boundary  Layers 

Estimates  have  been  made  of  the  propeller  blade  boundary-layer 
thickness  again  for  several  assumed  operating  conditions.  For  each 
of  these  conditions  it  is  estimated  that  less  than  5%  of  the  total  mass 
flux  of  fluid  passing  through  the  propeller  is  entrained  within  the 
propeller  blade  boundary  layers.  Thus,  even  if  the  turbulent  intensity 
within  the  blade  boundary  layers  were  appreciably  greater  than  that 
washing  off  the  body  and  through  the  propeller,  the  change  in  the  amount 
of  turbulent  energy  in  the  wake  would  be  small. 

If  the  blade  boundary  layers  were  to  separate  from  the  blade 
surfaces,  it  might  alter  the  above  conclusion.  It  should  be  recognized, 
however,  that  blade  boundary  layer  separation  would  probably  also 
result  in  markedly  reduced  propulsion  efficiency  since  the  thrust  (i.e., 
blade  lift)  decreases  and  the  torque  (i.e.,  blade  drag)  markedly  increases. 
Thus,  extensive  boundary  layer  separation  is  a pathological  case  corres- 
ponding to  operating  at  an  off-design  condition. 

Distortion  of  Turbulence  Passing  Through  a Marine  Propeller 

As  demonstrated  in  the  mean  flow  velocity  analysis  presented 
in  the  appendix,  the  primary  change  experienced  by  each  fluid  particle 
passing  through  the  propeller  is  that  it  is  suddenly  given  a swirl 
velocity  component.  The  changes  in  the  axial  and  radial  flow  components 
are  much  more  gradual.  This  imposition  of  a swirl  velocity  component  on 
the  turbulent  eddies  can  cause  their  orientation  to  change  and  they  may 
be  stretched  or  compressed.  Thus,  the  turbulence  within  the  aft-body 
boundary  layer  can  be  both  redistributed  into  other  components  and 
intensified  upon  passing  through  the  propeller.  Thus,  the  aft-body 
boundary-layer  turbulence  correlations  given  in  Section  2.1.2  could  be 
altered  by  virtue  of  the  propeller. 
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In  analyzing  this  change  it  Is  assumed  that  the  turbulence 
model  equations  presented  In  Section  3 (and  which  are  used  to  compute 
the  flow  aft  of  the  propeller)  remain  valid  even  in  the  vicinity  of  the 
propeller.  Thus,  these  equations  can  be  analytically  integrated  across 
the  plane  of  the  propeller  (i.e.,  from  z « -A  to  z - +A  , where  A is 
a very  small  distance  relative  to  the  propeller  or  wake  radius).  In 
doing  this,  it  is  assumed  that  each  of  the  turbulent  velocity  correlations 
such  as  Rzz  (Section  3)  could  experience  a step  change  from  its  value 
ahead  of  the  propeller  (say,  (RZZ)Q  )•  As  was  mentioned  above,  the 
mean  swirl  velocity  also  undergoes  a step  change,  AV,  being  zero  ahead 
of  the  propeller  and  finite  behind.  The  radial  distribution  of  the 
swirl  velocity  jump  is  computed  from  the  propeller  code  described  in 
Section  2.2.1. 

The  resulting  simultaneous  algebraic  equations  for  the  change 
in  each  turbulent  quantity  (i.e.,  R where  i , j ■ r , 0 , z ) 

and  the  turbulent  dissipation  rate  e can  then  be  solved.  Because  this 
is  rather  straightforward  (but  messy)  only  the  final  result  will  be 
given  as  follows  (see  Appendix  A. 3 for  details): 

Rrr  ‘ <Rrr)o  * °'582  <JT>2  F«>  . 

Ree  ' (Ree>o  “ l-°55  <r>2  F«>  - 

Rzz  ' <Rzz>o  - °-3M  <ir>2  F«>  • 

e - <£>„  [l  - 1.44  >2  *£>  ] , 

Rrz  * «rz>.  [2  - °-026  <iT>2]’ 

Rr0  ’ -°-236  <iT>  (Rrz>o  • 

Rez  ' f F(i>  ■ 

where  F(r)  - 0.327  (R  ) 0.018  (R.fl)  0.091  (R  ) , r - r/r  and 

zz  o Do  o rr  o p 

k2  m % J(RZZ)0  + ^R00^o  + (Rrr)0  j • Here,  r^  is  the  propeller  radius. 
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Note  that  changes  In  any  turbulence  quantity  are  solely  the  result  of 
the  sudden  Imposition  of  a swirl  velocity  to  the  stream  tube  passing 
through  the  propeller  plane. 

Most  of  these  changes  are  proportional  to  the  square  of  the 
normalized  swirl  velocity  with  the  exceptions  of  the  circumferential 
stresses  R^q  and  Rgz  which  are  proportional  to  Av  (the  constants 
shown  in  the  above  equations  result  from  the  turbulence  model  constants 
given  in  Section  3;  the  subscript  "o"  in  the  above  denotes  the  value 
of  that  quantity  immediately  upstream  of  the  propeller).  If  a typical 
value  of  (Av/w)  is  taken  as  20%  and  typical  values  of  the  boundary 
layer  correlations  used,  it  may  be  shown  that  the  changes  to  the  normal 
stress  components,  the  dissipation  rate  and  the  longitudinal  Reynolds 
stress  are  each  less  than  3%  of  their  values  upstream  of  the  propeller. 

The  swirling  Reynolds  stress  component  R^g  is  increased  from  zero 
upstream  to  about  5%  of  (Rrz)Q  » an<1  the  other  swirling  Reynolds  stress 
component  Rg^  is  changed  the  most,  going  from  zero  upstream  of  the 
propeller  to  about  90%  of  (Rrz)0  • Thus,  it  is  concluded  that  the 
propeller  has,  at  most,  a minimal  effect  on  the  turbulence  flowing 
through  it  except  for  causing  the  swirling  Reynolds  stresses  to  assume 
relatively  large  values. 

2.3  Section  Summary 

For  the  boundary-layer  mean  flow  near  the  aft-end  of  the  body 
a surprisingly  universal  correlation  of  Patel's  data  with  a flat-plate 
boundary  layer  has  been  described.  This  correlation  circumvents  the  problems 
associated  with  the  large  change  in  the  ratio  of  the  boundary-layer 
thickness  to  the  body  radius  and  with  the  strong  adverse  pressure  gradients 
present  near  the  trailing  edge.  A similar  correlation  of  the  boundary- 
layer  turbulence  properties  was  also  developed.  Thus,  given  the  body 
speed  and  size,  it  is  now  possible  to  make  estimates  of  the  aft-body 
boundary-layer  profiles  just  upstream  of  the  propeller. 

An  examination  of  the  flow  through  a marine  propeller  was  also 
conducted.  It  was  concluded  that  the  existing  Flow  Research  computer 
program  provides  a fairly  accurate  estimate  of  the  Induced  axial  and 
swirl  velocities.  An  alternate  approach,  presented  in  Appendix  A, 
provides  an  expression  for  the  induced  radial  velocity.  Adding  a routine 
to  compute  this  expression  to  the  existing  computer  program  provides 
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a complete  computer  model  for  the  mean  flow  in  the  propeller  plane.  The 
distortion  of  the  boundary- layer  turbulence  as  it  passes  through  the 
propeller  was  also  Investigated.  By  and  large,  the  changes  in  the 
turbulence  caused  by  the  propeller  appear  to  be  minimal  with  the  exception 
of  the  changes  in  the  swirling  Reynolds  stresses.  These  stresses,  because 
of  the  sudden  increase  in  swirl  velocity,  immediately  assume  values 
comparable  to  the  turbulent  shear  stress  within  the  upstream  boundary 
layer. 
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3.  Near-Wake  Region 

3.1  Mathematical  Model 

3.1.1  Mean-Flow  Equations 

In  the  near-wake  region,  we  assume  that  the  flow  is  statistically 
steady  and  axisymmetrlc,  and  that  the  effects  of  variations  in  density 
are  negligible.  We  perform  Reynolds'  decomposition  and  let  U , 

V and  W represent  mean  (time-averaged)  velocity  components  in  the 
r , 0 and  z directions  of  a cylindrical  polar  coordinate  system 
fixed  with  respect  to  the  body  (Figure  3.1)  with  its  origin  at  the  body 
tail.  The  z-axis  is  located  along  the  body  track.  We  denote  the 
fluctuating  velocities  by  u'  , v'  and  w'  and  correlations  of 
fluctuating  components,  that  is  u'u'  , u'v'  , etc.,  by  R^r  , 

Rre  , etc.  In  non-dimensional  form,  the  radial,  circumferential 
and  axial  momentum  equations  are 
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We  express  mass  conservation  as 


(3.1) 


(3.2) 


(3.3) 


(3. A) 


(3.5) 


Note  that  by  axisymmetrlc  flow  we  mean  that  circumferential  gradients 
of  mean  quantities  vanish. 
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The  mean  and  fluctuating  velocity  components  are  non- 
dimensionalized  by  a characteristic  mean  velocity,  and  lengths  by  a 
characteristic  radial  length  scale;  these  scales  are  th  free-stream 
velocity,  , and  the  body  radius,  R^  , respectively.  The  pressure, 

F , is  non-dimensionalized  by  pW^  , where  p is  the  local  fluid 
density,  after  subtraction  of  the  ambient  pressure  far  from  the  wake 
(so  that  P -+  0 as  r ■+■  °°).  The  flow  Reynolds  number.  Re  , is  defined 

as  WooVV- 

In  order  that  we  may  determine  the  pressure,  we  introduce  the 
Poisson  equation  for  P which  is  obtained  by  taking  the  divergence  of 
the  mean  momentum  equations.  This  results  in  an  equation  of  the  form 

2 2 

£ K 1 ^ P 

7T  + 7^  + 72  - F(u-  v-  "•  Rrr*  Ree'  «<=•>  «•<>> 

dr  d z 

Since  this  equation  is  used  in  place  of  the  mean-flow  continuity  equation 
(Equation  3.4),  the  formulation  and  form  of  the  function  F require 
careful  consideration.  We  must  ensure  that  P , as  numerically  pre- 
dicted by  this  equation,  be  such  that  the  corresponding  mean  velocity 
field  be  divergence  free. 

3.1.2  Turbulence  Model  Equations 

The  mean  flow  equations  require  the  specification  of  six 
turbulent  correlations,  Rrr  , R00  , R^  , RrQ  , Rfz  , and  R0z  , 

none  of  which,  in  general,  vanish  in  flows  with  swirl,  even  with  the 
assumption  of  axisymmetry.  To  obtain  these  correlations,  we  apply  the 
second-order  closure  turbulence  model  which  was  proposed  by  Hanjalic  and 
Launder  (1972),  and  extended  by  Launder,  Reece  and  Rodi  (1975).  Two 
factors  entered  our  decision  to  use  a second-order  closure  model. 

First,  and  perhaps  foremost,  was  the  lack  of  success  reported  by  a 
number  of  authors  (see,  for  example.  Pope  and  Whitelaw,  1976)  in  the 
application  of  simpler  turbulence  models,  such  as  the  several  k-k£ 
and  k-e  two-equation  models,  to  complex  turbulent  flows  which  are 
less  complex  than  flows  with  swirl.  First-order  closure  models  require 
constitutive  relationships  between  k and  kfc  (or  e ) and  the  six 
unknown  correlations;  typically,  the  Reynolds  stress  tensor,  T , is 
assumed  to  equal  yfidef  U , where  yg  is  an  isotropic  eddy  viscosity 
calculated  from  k and  kZ  (or  £ ).  Occasionally,  an  anisotropic 
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eddy  viscosity  is  generated,  although  on  purely  ad  hoc  grounds  (see, 
for  example,  Lilley,  1973).  For  complex  flows,  these  constitutive 
relationships  are  not  physically  well-founded. 

The  second  factor  derives  directly  from  the  motivation  for 
this  study.  We  are  attempting  to  develop  an  analytical/numerical  pro- 
cedure which  will  provide,  at  an  appropriate  station  downstream  of  the 
body,  physically  realistic  profiles  of  not  only  the  mean  flow  quantities 
(for  which  perhaps  the  use  of  a simpler  model  might  suffice),  but  also 
the  turbulence  quantities.  However,  two-equation  models  have  been 
developed,  in  general,  in  order  to  model  the  effects  of  turbulent  transport 
on  the  mean  floW — not  to  model  the  properties  of  the  turbulence  itself, 
and  thus,  we  do  not  expect  them  to  be  sufficient  for  our  purposes. 

In  terms  of  rectangular  cartesian  tensors  (and  assuming  the 
usual  summation  convention),  the  Hanjalic  and  Launder /Launder , Reece  and 
Rodi  equation  for  the  R. , turbulent  correlation  is 
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Where  the  terms  balancing  advectlon  of  the  correlation  represent,  in 
order:  production  resulting  from  the  interaction  of  the  fluctuating  and 
mean  flows,  isotropic  (high  Reynolds  number)  dissipation,  the  so-called 
"tendency  toward  isotropy"  component  of  the  pressure-strain  correlation 
(•t^j  ^),  the  mean  strain  component  of  the  pressure-strain  correlation, 
and  finally  turbulent  diffusion.  The  mean-strain  component  of  the 
pressure-strain  correlation  (as  developed  by  Launder,  et  al.)  is  defined 
as 


ij,2 


(C2  + 8) 

IT 


2 (30Cj  - 2)  / 3U,  3U 

(Pij  * 3 Pk6ij)  55  k \ 3^i  ’ 3Xj 


(8C2  - 2) 
li 


<Du  - f W 


> 


I 


» 


I 


I 


I 


( 


I 


( 


c 


Flow  Research  Report  No.  69 
December  1978 


where 

P 


ij 


■(' 


3U 


‘1)1  3x. 


+ R 


3U 

V l 3^ 


i\  / 3Ui  3V 

-j  • and  Du  -■  (Ru^  + V^ 


and  represents  the  rate  of  turbulent  energy  production,  JjP^. 

The  two  underlined  turbulent  diffusion  terms  are  neglected  in  the  often- 
used  diffusion  expression  of  Daly  and  Harlow  (1970).  The  quantity  k 
is  the  turbulent  kinetic  energy  which  is  equal  to  IjR^  , and  £ is  the 
turbulent  energy  "dissipation  rate",  determined  from  the  model  equation, 


De 

Dt 


C — R C — + c — / — R.  \ 

k if.  9x^  k t 3x^  \ £ 3x^  ) 


(3.8) 


Here,  the  terms  balancing  advection  of  dissipation  rate  represent  its 
production,  dissipation  and  turbulent  diffusion. 

The  correlation  and  dissipation  rate  equations  are  presented 
in  terms  of  the  cylindrical  coordinate  system  of  this  study  in  Appendix 
B.  The  model  constants  are  also  defined  in  that  appendix.  As  a check 
of  the  equations,  we  note  that  for  non-swirling  flows,  V - 0 and, 
as  a result,  R^  **  Rgz  ■ 0 ; with  the  simplified  Daly-Harlow  expression 
for  diffusion,  the  model  equations  in  cylindrical  coordinates  then 
reduce  identically  to  those  used  by  Pope  and  Whitelaw  (1976). 

As  an  inspection  of  the  equations  in  Appendix  B will  verify, 
the  model  equations  are  exceedingly  lengthy.  This  is  due,  in  part, 
to  the  modeling  assumed  by  Hanjalic  and  Launder  for  the  turbulent 
diffusion  terms.  A rather  demanding  constraint  on  the  time  period 
of  this  research  effort  suggested  the  application  of  a simplifying 
hypothesis — we  assume  that  derivatives  with  respect  to  z in  the 
turbulent  diffusion  terms  may  be  neglected.  This  is  a boundary- layer- 
like  assumption  which  we  expect  will  have  a minimal  effect  on  the 
prediction  of  the  non-recirculating  wake  flows  of  interest  here. 

The  inclusion  of  these  terms  should  be  carefully  considered,  of  course, 
if  a model  such  as  this  were  to  be  applied  to  such  recirculating  flows 
as  highly  swirling  jets.  With  this  simplifying  assumption,  the  model 
equations  become  first-order  in  z ; they  are  second-order  in  r only. 
The  simplified  model  diffusion  terms  are  included  in  Appendix  B. 

3.1.3  Boundary  Conditions 


( 


The  system  of  equations  presented  in  the  previous  section  is 
elliptic  and  boundary  conditions  are  required  on  all  boundaries  of  the 
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solution  domain  (see  Figure  3.1).  In  what  follows  (and  in  the  computer 

code  itself),  boundaries  and  orientations  are  given  in  terms  of  the 

directions  North,  South,  East  and  West.  The  North  boundary  is  located 

at  radial  position,  r = R , and  the  South  boundary  is  at  the  axis  of 

max 

symmetry  (r  = 0).  The  West  Boundary,  at  z ■ , is  a plane  transverse 

to  the  axis,  located  just  aft  of  the  body  tail,  and  the  East  boundary  is 

a transverse  plane,  at  some  distance,  Z , downstream.  Thus  the 

domain  in  which  the  solution  is  to  be  obtained  is  defined  by  0 < r < R , 

— — max 

< z < Z 
i — — max 

At  the  axis  of  symmetry  (South  boundary),  continuity  requires 

that  U = 0 . Then  from  symmetry  conditions  along  with  an  examination  of 

the  mean  flow  equations  (keeping  in  mind  that  shear  stress  must  vanish 

at  the  axis),  we  find  that  V , 3W/3r  , 3P/3r  , 3R  /3r  , RQQ/3r  , 3R  /3r  , 

rr  ot)  zz 

3e/3r  , RfQ  , R^z  , and  Rgz  must  all  vanish.  Under  these  conditions, 
the  turbulence  model,  as  it  should,  will  predict  R^  » Rgg  at  the  axis 
(Pope,  1976). 

We  assume  that  the  radial  (North)  boundary  is  located  sufficiently 
far  removed  from  the  axis  of  symmetry  (and  body  track)  and  that  all 
turbulence  quantities, = as  well  as  the  propeller  induced  swirl  velocity, 
vanish.  We  assume  also  that  the  axial  velocity  is  uniform  and  equal  to 
its  free-stream  value;  then,  since  3W/3z  ■ 0 , we  allow  a radial  flux 
into  the  domain  with  the  condition  that  3(Ur)/3r  • 0 , which  is  obtained 
from  the  continuity  equation.  We  obtain  the  radial  pressure  gradient, 

3P/3r  , at  Rmax  from  the  mean-flow  radial  momentum  equation,  and  apply 
it  as  a boundary  condition  on  P . 

At  the  West  (inflow)  boundary,  we  assume  that  the  mean  velocity, 
turbulent  correlations,  and  dissipation  rate  are  all  specified.  Those 
quantities  which  are  not  experimentally  or  theoretically  known  must  be 
estimated.  We  compute  the  axial  pressure  gradient,  3P/dz  , from  the 
axial  mean-flow  momentum  equations,  and  apply  it  as  the  boundary  condition 
on  P . 

As  a result  of  the  simplifying  approximation  discussed  in  the 
previous  section,  no  conditions  are  required  for  the  turbulence  quantities 
at  the  East  boundary;  however,  we  must  supply  conditions  on  the  mean 
velocity  components  and  the  pressure.  The  conditions  applied  in  this 
study  are  given  below: 
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(3.10) 


(3.11) 


(3.12) 


This  system  of  equation  is  valid  for  slender,  swirling  shear 
flows  in  which  radial  gradients  are  in  general  much  larger  than  axial 
gradients,  and  it  is  sometimes  known  as  the  "quasi-cylindrical  approxi- 
mation" (Hall,  1966).  It  is  simply  a boundary-layer  approximation  in 
which  P , Instead  of  being  uniform  across  the  flow,  varies  to  balance 
centrifugal  acceleration  and  turbulent  transport.  This  system  has 
been  applied  by  many  investigators  (e.g.,  Hall,  1965;  Mager,  1972) 
although  it  has  not  previously  been  applied  as  a boundary  condition. 

By  using  it  as  a boundary  condition,  we  avoid  the  necessity  of  applying 
much  cruder  approximations  (such  as  the  often-used  requirement  that 
axial  derivatives  of  flow  quantities  vanish) , and  our  solutions  will 
automatically  match  with  the  nonstratif led  far-wake  solution  which 
satisfies  these  equations. 
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3.1.4  Coordinate  Transformation 

In  order  that  we  may  choose  the  boundary  locations  R and 

max 

^max  sufficiently  large  that  the  above  boundary  conditions  may  be 
accurately  applied,  while  at  the  same  time  maintaining  sufficient  radial 
mesh  resolution  within  the  wake  (r  < 1) , and  axial  mesh  resolution, 
particularly  in  the  upstream  portion  of  the  wake  where  axial  gradients 
are  expected  to  be  greatest,  we  perform  logarithmic  coordinate  transforma- 
tions. We  introduce  the  independent  variables  y and  x where 

# » 

y m~  in(l  + r/b  ) , and  x - ~ £n(l  + (z  - *±)/b  ) . 

y y x 

These  have  been  applied  in  conjunction  with  numerical  solutions 
to  a number  of  flow  problems  (e.g.,  Pao  and  Daugherty,  1969,  and  Grabowski 
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and  Berger,  1976),  and  are  used  here  to  map  the  region  0 £ r £ R x » 

Z,  < z < Z , in  which  we  wish  to  obtain  a solution,  onto  the  region 
i — — max 

0 £ x £ 1 in  transformed  space.  The  values  of  the  four 

transformation  parameters  a , a , b , and  b are  determined  from 

y x y x 

the  specified  values  of  R and  Z , and  the  desired  number  and 
r max  max 

distribution  of  grid  points  in  each  of  the  coordinate  directions. 

Details  are  presented  in  Appendix  C.  Derivatives  of  dependent  variables 
transform  into  the  new  coordinate  system  as  follows: 
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f (y)  - exp(-ayy)/(ayby)  , 


g(y)  * ~ exp(-2ayy)/(ayby)  , 


h(y)  - 1/r (y)  , 


s(x)  - exp(-axx)/(axbx) 


and 


t(x)  - -exp(-2axx)/(axbx)  . 


In  terms  of  the  independent  variables  x and  y , the  trans- 
formed mean-flow  momentum  equations  are 
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The  equations  for  the  turbulent  correlations,  the  dissipation  rate, 
and  the  boundary  conditions  transform  similarly  and  are  not  shown  here. 

3. 2 Numerical  Procedure 

3.2.1  Introduction  of  T 

In  order  to  obtain  a solution  of  Equations  (3.13)  through  (3.15), 

along  with  the  transformed  turbulence  model  equations  and  the  boundary 

conditions,  we  introduce  a time-like  independent  variable,  T , 

and  add  the  "time  derivative"  terms  3U/3t  , 3V/3t  , 3W/3t  , 3e/3x  , 

3Rrr./3t  , 3Rqq/3t  , etc.,  to  Equations  (3.13)  through  (3.15)  and  the  model 

equations  for  e , R , RQQ  , etc.,  respectively.  We  also  modify 

rr  uu 

the  East  boundary  conditions  with  the  addition  of  3V/3x  and  3W/3x  to 
the  transformed  version  of  Equations  (3.10)  and  (3.11),  respectively. 

The  East  boundary  condition  thus  becomes  what  might  be  termed  an 
"unsteady  quasi-cylindrical"  approximation.*  The  introduction  of  a 
time-like  dependence  to  the  system  has  two  important  advantages.  First, 
it  allows  the  system  to  be  solved  to  the  steady  state,  in  which  the 
"time  derivatives"  vanish,  by  any  one  of  the  many  numerical  schemes 
developed  for  initial  value  problems.  The  steady  state  solution  is 
obviously  the  solution  we  desire.  The  conversion  of  a boundary-value 
problem  to  an  initial/boundary-value  problem  in  order  to  facilitate 
numerical  solution  Is  a standard  procedure.  We  also  gain  a second 
related  benefit  from  the  Introduction  of  a x-dependence  In  that  the 
application  of  the  East  boundary  condition  is  significantly  facilitated. 

The  equations  at  the  East  boundary  and  those  in  the  Interior  of  the 
domain  are  marched  forward  In  X In  unison  in  straightforward  fashion. 

We  thus  have  the  not  unusual  situation  of  an  Inltlal/boundary-value 
problem  with  a boundary  condition  which  varies  In  "time." 


Note  that  3U/3x  Is  not  added  to  Equation  (3.9).  For  the  purposes  of 
this  study  we  may  assume  that  3U/3x«3P/3r  at  the  East  boundary. 
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3.2.2  The  Pressure  Equation 

As  noted  in  an  earlier  section,  we  obtain  a Poisson  equation 
for  P by  taking  the  divergence  of  the  momentum  equations.  In  terms  of 
the  transformed  variables,  we  write  the  equation  as 


f2  M + (g  + fh)  + s2  + t - F(U,  V,  W,  R , etc.)  , (3.16) 

3y2  3y  3x2  9x  rr 


where 


, 3C  7 3C 

F - r -g*  + (g  + fh)Cy  + + tCx  + 3Z?/3x  , 
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-fU  - sW  + hV2  + 
3y  3x 
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-s  I -fU  37  - sW  ^ 


3R  3R 

zz  , , rz  , 

+ 8 —5 — + f — ^ — + hR 
3x  3y  rz 


and 


n it  ^ . -liT\  j. 

D m hf  37  01  U)  + 8 3^  • 


The  quantity  D Is  the  divergence  of  the  velocity  field  in  terms  of 
the  transformed  variables.  Although  continuity  (that  is,  the  condition 
D should  equal  zero)  has  been  applied  to  eliminate  the  viscous  terms, 
the  "time"  derivative  of  D must  be  retained  in  order  to  assure  stable 
solution.  The  necessity  of  retaining  3Z?/3t  when  the  Poisson  equation  for 
pressure  is  solved  numerically  was  first  noted  by  Harlow  and  Welch  (1965). 
The  motivation  for  the  retention  cf  this  term  is  that  in  the  finite- 
difference  computation  truncation  and  round-off  error  which  exist  at  any 
time  level,  n-^  (we  use  half-time  steps  in  the  finite-difference 
procedure),  result  in  a non-zero  divergence  which,  unless  properly 
treated,  grows  and  gives  rise  to  a computational  instability.  However, 
since  at  time-level  n , D ought  to  equal  zero,  we  approximate 
32?/3t  at  any  grid  point  by  finite-difference  at  time  level  n as  (0- 
Z?”_,5)/(1/2At)  . This  stratagem,  which  prevents  the  growth  of  D , was 
originally  suggested  by  Harlow  and  Welch. 

There  are  many  possible  mathematically  consistent  forms  in 
which  the  remaining  terms  in  F may  be  expressed.  The  formulation  used 
here  was  arrived  at  after  some  experimentation  and  consideration.  In 
essence,  each  term  in  F corresponds  to  a term  on  the  left-hand  side  of 
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the  equation  and,  in  fact,  is  arrived  at  by  using  the  appropriate 
equation  to  evaluate  the  derivatives  of  P . For  example,  the  term  tC^ 
is  obtained  by  using  the  x-momentum  equation  to  express  3P/3x  in 
terms  of  velocity  components  and  turbulent  correlations.  This  expression 
is  included  in  F to  correspond  to  the  t3P/3x  term  on  the  left.  The 
other  terms  in  F are  obtained  similarly.  We  note  that  we  have  combined 
all  of  the  "time  derivatives"  into  3Z?/3x  and,  since  the  Reynolds 
number  is  assumed  large  so  that  the  error  arising  from  the  finite- 
difference  representation  of  the  viscous  terms  is  negligible,  we  have 
used  continuity  (that  is,  the  assumption  that  D is  identically  zero) 
to  combine  and  eliminate  them. 


3.2.3 


Finite-Difference  Formulation  and  Numerical  Procedure 


To  generate  a numerical  procedure  for  the  solution  of  Equations 

(3.13),  (3.14),  (3.15),  (3.16),  and  the  seven  transformed  turbulence 

model  equations,  we  first  overlay  the  region  0 £ x ^ 1 , 0 .£  y .1  *5  , 

in  which  the  solution  is  to  be  obtained,  with  a finite-difference  grid 

system  of  N points  in  the  x-direction  and  M points  in  the  y- 

direction.  The  location  (x,y)  of  any  grid  point  is  given  by  the 

indices  i,j  as  x - (i  - l)Ax,  y « (j  - l)Ay  , where  Ax  and  Ay 

are  uniform  (but  not  equal)  mesh  widths  in  the  x and  y directions, 

respectively  (Figure  3.2).  The  "unsteady"  solution  functions  U(x,y,t)  , 

V(x,y,T)  , W(x,y,T)  , Rrr(x,y,x)  , etc.,  are  represented  by  the  discrete 

functions  u"  , v”  , W?  , Rn  , etc.,  which  are  defined  at  each 
ij  ij  ij  rr;Lj 

grid  point  i,j  . The  superscript  n refers  to  the  number  of  discrete 
time  steps  of  length  Ax  equivalent  to  x , i.e.  Xs  nAx  . The  pressure 
P(x,y,x)  is  represented  by  the  discrete  function  P^^  which  is 

defined  on  a second  grid  system  staggered  with  respect  to  the  first  by 
half  of  a mesh  width  in  both  coordinate  directions. 

The  finite-difference  formulation  of  the  transformed  mean  flow 
equations  (Equations  3.13  through  3.15)  is  based  on  the  well-known  ADI 
(Alternating-Direction- Implicit)  procedure  of  Peaceman  and  Rachford 
(1955),  and  the  finite-difference  equations  are  presented  in  detail 
in  Appendix  D.  We  approximate  all  derivatives  in  the  viscous  diffusion, 
Reynolds  stress,  pressure  gradient  and  radial  convection  terms  with  the 
usual  centered  second-order  finite-difference  approximations,  while 


A 


c 


Flow  Research  Report  No.  69 
December  1978 


-29- 


we  approximate  axial  advectlon  terms  such  as  WdU/dx  with  a weighted 
centered-upwlnd  formula  of  the  form 

[ui  ♦ i.jO'u  - + ' ui  - + «iKjl)],2te  • 

The  values  permitted  for  the  variable  weighting  factor  are  bounded 

by  zero  and  unity  in  which  cases  the  formula  yields,  respectively, 
second-order  centered  and  first-order  upwind  differencing. 

The  ADI  procedure  Introduces  two  halftime-steps  of  length 

4(At)  , during  which  y and  x derivatives  of  the  mean-flow  velocities 

U , V , and  W , including  the  derivatives  in  the  convection  terms,  are 

alternately  differenced  implicitly  (in  time).  As  a result  of  this 

n 1 1,-  n+J< 

procedure,  three  tridiagonal  matrix  equations  for  * , 

W^*4  , are  defined  at  each  value  of  i for  2 £ j £ M during  the  first 
halftime-step.  When  modified  to  incorporate  the  appropriate  boundary 
conditions,  these  equations  yield  the  solution  at  the  grid  points  i , 

1 £ 3 £ M • During  the  second  half-step,  three  tridiagonal  equations 
are  defined  at  each  value  of  j , and  these,  after  incorporation  of  the 
appropriate  boundary  conditions,  yield  the  solution  , and 

at  the  grid  points  j , 2 £ i < N (recall  that  the  solution  is 
specified  at  i » 1).  The  tridiagonal  matrix  equations  are  solved  using 
the  simple  Thomas  variant  (Roache,  1972)  of  Gaussian  elimination  with  a 
modification  which  allows  the  incorporation  of  Neumann  boundary  conditions 
to  second-order. 

Since  the  seven  equations  which  determine  the  turbulence 
quantities  contain  second-order  derivatives  in  y only,  implicit  dif- 
ferencing is  applied  only  in  the  y direction.  Centered  spatial 
differences  are  used  to  approximate  all  derivatives  except  those  appear- 
ing in  axial  advectlon  terms.  As  described  above,  we  approximate  these  with 
weighted  centered-upwind  differences.*  During  every  half  time-step, 
seven  tridiagonal  matrix  equations,  into  which  the  North,  South  and  West 
boundary  conditions  are  incorporated,  are  inverted  for  each  value  of  1 
(except  i * 1)  to  yield  the  turbulence  quantities  at  the  grid  points  i , 

1 < J < H . 


The  finite-difference  representation  of  the  equation  for  R is 
included  In  Appendix  D. 
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The  discretized  pressure,  P , is  defined  on  the  staggered  grid 
system  denoted  by  i+ij,  j+ij  in  order  to  simplify  the  application  of  the 
Neumann  pressure  boundary  conditions  at  the  North,  South  and  West 
boundaries  and  to  Increase  the  accuracy  of  the  scheme.  Centered  second- 
order  differences  are  used  to  approximate  all  of  the  derivatives  of  P 
in  Equation  3.16.  Centered  differences  are  also  used  to  approximate  all 
of  the  derivatives  appearing  in  the  expression  F except  for  the  terms 
sW3U/3x  and  sW3W/3x  which  appear,  respectively,  in  C and  C 

y * 

These  advection  terms  are  approximated  with  the  weighted  centered-upwind 
formula  described  above.  The  finite-difference  approximation  to  the 
pressure  equation  is  included  in  Appendix  D. 

To  obtain  the  solution  of  the  discretized  pressure  equation  at 
each  halftime-step,  we  employ  the  direct  (that  is,  noniterative)  technique 
of  Swartztrauber  (1974).  The  method  is  based  on  the  cyclic  reduction 
technique  and  is  fully  described  in  Swartztrauber' s paper.  It  was 
developed  for  the  numerical  solution  of  general  second-order  separable 
elliptic  partial  differential  equations  in  which  the  coefficients  or 
differential  terms  vary  with  the  particular  variable  of  differentiation.* 
The  transformed  pressure  equation  derived  in  this  study  is  an  excellent 
example  of  such  an  equation,  and  application  of  the  method  is  straight- 
forward. We  were  fortunate  to  obtain  a copy  of  the  computer  code  BLKTRI, 
which  is  based  on  the  method,  from  the  National  Center  for  Atmospheric 
Research.  The  routine  BLKTRI  is  incorporated  into  the  ICWAKE  code. 

Each  solution  of  the  (N  - 1) (M  - 1)  interior  values  of  P requires 
only  (M  - 1)(N  - lHogj  (M  - 1)  operations  excluding  the  operations 
required  to  compute  F . The  procedure  is  thus  quite  efficient. 

The  transformed  counterparts  of  Equations  (3.10)  and  (3.11) 
which  apply  at  the  East  boundary  are  approximated  with  central  dif- 
ferences for  all  derivatives  except  those  in  the  axial  convection  terms 
which  are  approximated  with  first-order  backward  differences.  During 
the  alternating  halftime-steps,  the  y and  x derivatives  are  treated 
Implicitly.  Two  tridiagonal  matrix  equations  result  during  the  first 
half  step  for  * T*le8e  may  be  inverted,  and  the  solution 

* 

We  note  that  the  procedure  may  also  be  applied  to  the  matrix  problem 
generated  by  the  discretization  of  a separable  elliptic  equation  over 
a non-uniform  finite-difference  grid  with  mesh  spacing  which  varies 
independently  in  the  two  coordinate  directions.  We  also  note  that  by 
separable  we  mean  that  the  differential  operator  Is  separable.  The 
inhomogeneous  terms  need  not  be  separable. 
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along  with  the  appropriate  turbulence  quantities  are  used  to  Integrate 


Equation  (3.9)  for  pH?1V  and  Equation  (3.12)  for  . These  values  of 

n+k  nj  n+i 

?Nj+4  are  use<*  as  boundary  conditions  for  the  calculation  of  P In 

the  interior  of  the  domain.  During  the  second  half-time  step,  Equations 

(3.10),  (3.11),  and  (3.12)  are  incorporated  directly  into  the  matrix 

equations  for  each  j . After  these  equations  are  inverted.  Equation 

(3.9)  is  used  to  obtain  P^V  which  is  used  as  a boundary  condition  for 

the  calculation  of  pn+^  . 

As  implied  in  the  above  discussion,  a calculation  through  a 

complete  time  step  proceeds  as  follows.  Based  on  the  known  values  of 

n*f^s  n+^ 

all  the  flow  variables  at  time  step  n , the  values  , 

and  are  first  obtained  for  all  i,j  from  the  momentum  equations 

with  y derivatives  treated  implicitly.  With  these  values,  the  seven 

turbulence  quantities  are  then  obtained  at  n + *s  . To  complete  the 

n-fk 

half-step,  the  pressure  P is  calculated  using  the  function  F 
evaluated  at  n+lj  with  the  known  values  of  the  other  variables.  The 
second  time-step  essentially  repeats  the  first  and  yields  the  solution 
at  time  level  n + 1 . The  difference  between  the  first  and  second  time 
half-steps  is  in  the  treatment  of  the  mean-flow  momentum  equations  only. 

It  should  be  obvious  that  this  finite-difference  procedure  can 
be  used  for  the  calculation  of  laminar  flows.  For  such  a calculation, 
the  computation  of  the  turbulent  correlations  is  simply  not  performed. 
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4.  Sample  Calculations 

In  this  section,  we  present  results  from  four  near-wake 
calculations.  The  first  calculation  is  initialized  by  and  compared  with 
the  laboratory  data  of  Schetz  (reported  by  Swanson,  et  al.  (1974)  and 
Chieng,  et  al.  (1974)).  The  second  calculation  is  initialized  by  and 
compared  with  the  data  of  Chevray  (1968),  and  is  also  compared  with  a 
calculation  of  Pope  and  Whitelaw  (1976).  The  third  and  fourth  calcula- 
tions are  initialized  with  profiles  obtained  from  the  boundary  layer/ 

A 

propeller  model  described  in  Section  2.  The  third  is  for  Re  - 1.5  * 10 
and  is  compared  with  the  data  of  Lin,  et  al.  (1974)  at  z/R^  “ 12  . The 
fourth  is  for  Re  * 1.5  x 10^  . Comparative  plots  of  wake  width  versus 
z for  the  third  and  fourth  calculations  are  shown. 

4.1  Calculation  Initialized  By  the  Data  of  Schetz 

The  data  of  Schetz  were  obtained  in  the  wake  behind  a blunt- 

nosed,  parallel-sided,  sharp-stemed,  axisymmetric  body  with  a maximum 

diameter  of  6 inches  and  a length  of  60  inches.  The  flow  Reynolds 

number,  based  on  free-stream  velocity  and  body  radius,  was  3.1  * 10^  . 

The  body  was  fitted  with  a three-bladed  6-inch-diameter  aft-mounted 

propeller  and  the  experiments  were  conducted  under  approximately  drag- 

free  conditions.  Hot  wire  data  at  axial  positions  4,  10,  20,  40  and  80 

body  radii  downstream  of  the  body  tail  were  obtained.  These  data 

include  the  turbulence  correlations  R , RQD  , R , R.  , and 

zz  oo  rr  oz 

Rrz  , as  well  as  W and  V , the  axial  and  swirl  mean  velocity  components. 

In  the  ICWAKE  calculation  presented  here,  the  West  boundary 
was  set  at  the  first  data  station,  4 body  radii  behind  the  body,  and  the 
East  boundary  was  set  at  the  third  data  station,  20  body  radii  down- 
stream. The  North  boundary  was  2 body  radii  from  the  axis  of  symmetry. 

The  code  was  initialized  at  the  West  boundary  by  the  seven 
measured  profiles  and  the  following  estimates  for  the  unmeasured  quantities 
U , R^g  , and  e . The  mean  radial  velocity  component,  U , was  set  to 
zero;  the  turbulence  correlation  R^g  was  obtained  from  the  eddy 
viscosity  approximation 
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where 


V -R  |S 

e rz  9r 


and  the  dissipation  rate  £ was  set  equal  to  K k /A  , where 

2 1 e 
k - -r-(R  + Rqq  + R ),  K is  a constant,  and  A is  the  integral 

/ rr  ot)  zz  £ 

length  scale  of  high  turbulent  Reynold's  number  theory.  The  values  of 
and  A were  chosen  to  be  0.53  and  0.2  (i.e.,  0.2  body  radii),  as 
suggested  from  the  experimental  data  obtained  by  Gran  (1973). 

The  grid  system  for  the  calculation  consisted  of  40  points  in 
the  axial  direction  with  10  of  the  points  between  z/R^  * 4 and 
z/Rj^  = 5 , and  32  points  in  the  radial  direction  with  16  of  the  points 
between  r/R^  = 0 and  r/R^  = 0.8  . Convergence  was  assumed  when  the 
rate  of  change  of  the  RMS  mean  velocity  divergence,  calculated  over  the 

-9 

computational  grid,  became  less  than  10  . The  time-step  used  was  0.1; 

about  six  minutes  of  CDC  7600  central  processor  time  was  required  to  do 
the  calculation. 

Comparisions  of  measured  and  computed  cross-wake  profiles  at 

z/R^  *=  10  are  presented  in  Figures  4.1  through  4.6.  Figure  4.1  displays 

the  axial  and  swirl  mean  velocity  components.  The  predictions  are  in 

excellent  agreement  with  the  experimental  data.  Figures  4.2,  4.3,  and 

4.4  present  the  three  normal  correlations  R , RQQ  , and  R , 

zz  ow  rr 

respectively.  The  predicted  values  of  R are  in  reasonable  agreement 

ZZ 

with  the  experimental  data  although  the  secondary  peak  in  the  data  at 
about  r/R^  “ 1 * which  is  presumably  due  to  the  vortex  system  shed  from 
the  propeller  tips,  does  not  appear  in  the  prediction.  The  peaks  in 
Rqq  and  Rrr  measured  near  the  wake  edge  also  do  not  appear  in  the 
predictions  of  these  quantities,  and  in  the  case  of  R^^  , the  difference 
between  the  predicted  and  measured  values  is  quite  large,  although  the 
value  of  Rrr  predicted  at  the  axis  is  in  good  agreement  with  the  data. 
Figures  4.5  and  4.6  display  the  radial  and  tangential  shear  stresses, 

Rfz  and  Rqz  , across  the  wake.  The  agreement  between  the  predictions 
and  the  data  for  Rfz  is  quite  good  throughout  the  wake,  and  for  Rqz 
is  reasonable  in  the  interior  of  the  wake.  Near  the  edge  of  the  wake, 
however,  the  large  negative  measured  values  of  R0  are  not  predicted. 


are  not  oredicted. 


4.2  Calculation  Initialized  By  the  Data  of  Chevray 

The  data  of  Chevray  were  obtained  in  the  wake  behind  a five 
foot  long,  slx-to-one  prolate  spheroid.  The  flow  Reynolds  number,  based 
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on  free-stream  velocity  and  body  radius,  was  2.3  * 10J  , and  the  drag 

coefficient  was  0.06.  Hot  wire  data  at  axial  positions  of  -0.5,  0,  0.5, 

1.,  2.,  4.,  6.,  12.,  18.,  24.,  30.,  and  36  body  radii  downstream  of  the 

body  tail  were  obtained.  These  data  include  the  turbulence  correlations 

R , R00  , R , and  R (R„  ■ RQ  = 0)  , and  the  mean  axial  and 

zz  ou  rr  rz  oz  or 

radial  velocity  components  W and  U (V  = 0)  . 

In  the  ICWAKE  calculations  presented  here,  the  West  boundary 
was  set  at  z/R^  = 0.5  and  the  East  boundary  was  set  at  z /R^  * 20  . 

The  North  boundary  was  4 body  radii  from  the  axis  of  symmetry. 

The  calculation  was  initialized  at  the  West  boundary  by  the 
six  measured  profiles,  the  knowledge  that  R^z  » R^r  * V * 0 , and  the 
following  assumptions  about  the  dissipation  rate  e . From  the  North 
boundary  to  r/R^  = r*  » where  the  production  of  turbulent  energy,  P , 
assumes  its  maximum  value,  Pq  , e was  set  equal  to  P . From  r/R^  = r* 
to  the  axis  of  symmetry,  e was  set  equal  to  Pq  . 

The  grid  system  for  the  calculation  consisted  of  40  points  in 
the  axial  direction  with  10  of  the  points  between  z/R^  - 0.5  and 
z/R^  * 1.5  , and  32  points  in  the  radial  direction  with  18  of  the 
points  between  r/R^  “ 0 and  r/R^  ” 1 • The  convergence  criterion  was 
the  same  as  in  the  Schetz  case.  A time-step  of  0.05  was  required  for 
convergence,  and  about  10.5  minutes  of  CDC  7600  central  processor  time 
was  needed  to  perform  the  computation. 

Comparisons  of  measured  and  computed  cross-wake  profiles  of 
W , R^^  , and  R^  at  z/R^  * 6 , 12  , and  18  are  presented  in 
Figures  4.7  through  4.9.  For  W and  Rfz  , the  ICWAKE  results  are  also 
compared  to  the  published  computational  results  of  Pope  and  Whltelaw 
(1976).  Figures  4.7a,  b,  c present  the  comparisons  for  the  axial  mean 
velocity,  W . The  ICWAKE  calculation  is  in  excellent  agreement  with 
the  Pope  and  Whltelaw  calculation,  and  both  calculations  give  good 
predictions  of  Chevray's  data.  The  predictions  are  slightly  less  reliable 
at  the  axis  of  symmetry  than  near  the  wake  boundary,  but  this  discrepancy 
is  eye-catching  only  at  the  downstream  station,  z/Rb  » 6 . Figures 
4.8a,  b,  c display  the  results  for  the  radial  shear  stress.  Again  the 
ICWAKE  results  agree  very  closely  with  those  of  Pope  and  Whltelaw. 

Although  both  calculations  overestimate  the  peak  magnitudes  of  the  shear 
stress,  they  give  reasonable  predictions  of  the  data.  Lastly,  Figures 
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A. 9a,  b,  c show  comparisons  between  ICWAKE  results  and  Chevray's  data 
for  the  axial  turbulence  correlation.  Results  for  this  variable  were 
not  presented  by  Pope  and  Whitelaw.  The  ICWAKE  results  overestimate  the 

I 

data  values,  but  show  the  basic  trends.  The  largest  overpredictions  are 
at  the  axis  of  symmetry,  and  the  best  predictions  are  near  the  wake 
edge. 

A. 3 Calculation  Initialized  By  the  Boundary-Layer /Propeller  Model. 

Re  = 1.5  x IQ4 

The  previously  described  boundary-layer/propeller  model  was 
used  to  simulate  a Flow  Research  laboratory  experiment  described  by 
Lin,  et  al.  (197A).  In  this  experiment,  data  were  obtained  in  the  wake 
behind  a body/propeller  configuration  similar  to  that  of  Schetz.  The 
maximum  body  diameter  was  10.16  cm,  the  body  length  was  119  cm,  the 
propeller  diameter  was  A. 3 cm,  and  the  Reynold's  number.  Re  = 1.5  * 10^  . 

Hot  film  data  at  axial  positions  12,  30,  A8,  and  1AA  body  radii  down- 
stream of  the  body  tail  were  obtained.  These  data  included  the  turbulence 

correlations  R , R , and  R , and  the  axial  mean  velocity 
zz  4 r rz 

component  W . 

In  the  ICWAKE  calculation,  the  West  boundary  was  set  ac 

z/R^  » 0 where  all  initial  conditions  were  specified  by  the  boundary- 

layer  /propeller  model.  The  propeller  parameters  required  by  the  model 

were  obtained  from  measurements  of  the  propeller  on  the  laboratory 

vehicle.  The  East  boundary  was  set  at  20  body  radii  downstream  and  the 

North  boundary  was  set  at  2 body  radii  from  the  axis  of  symmetry. 

Because  of  numerical  problems  which  were  encountered,  special 

treatment  of  the  initial  turbulence  intensity  and  dissipation  rate 

profiles  was  required.  The  predicted  profiles  of  these  variables  at 

the  initial  station,  z/R^  ■ 0 , were  such  that  the  variable  values 

dropped  to  zero  roughly  halfway  between  the  axis  of  symmetry  and  the 

edge  of  the  wake.  As  experience  showed,  the  code  would  fall  catastrophically 

unless  these  profiles  were  non-zero  throughout  the  wake.  So  from  the 

point  at  which  these  variables  became  sufficiently  small  to  be  dynamically 

insignificant,  say  less  than  10  to  the  edge  of  the  wake  where  they 

were  constrained  by  the  North  boundary  conditions  to  be  zero,  a linear 

approximation  was  used.  In  addition,  it  was  found  necessary  to  assign 

f— 2 5 ~Tl 

an  upper  bound  to  eW^Rzz  + rqq  + Rrr)  » In  the  outer  wake  region.  This 
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bound  was  chosen  from  examination  of  previous  successful  calculations  to 
be  10  . 

The  grid  system  for  the  calculation  consisted  of  60  points  in 
the  axial  direction  with  20  of  the  points  between  z/R^  ■ 0 and 
z/R^  ■ 1 , and  32  points  in  the  radial  direction  with  24  of  the  points 
between  r/R^  “ 0 and  r/R^  * 1.2  . The  convergence  criterion  was 
again  as  in  the  Schetz  case.  A time-step  of  0.05  was  used,  and  about 
eleven  minutes  of  CDC  7600  central  processor  time  was  required  to  do  the 
calculations. 

Comparisons  of  measured  and  computed  cross-wake  profiles  at 
r/Rb  = 12  are  presented  in  Figures  4.10  and  4.11.  The  quantities 
displayed  are  W and  Rzz  , respectively.  The  agreement  between 
predictions  and  measurements  is  quite  reasonable  in  both  cases  except  for 
the  axial  mean  velocity  component  near  the  axis  of  symmetry. 

Figure  4.12  shows  a plot  of  the  width  of  the  simulated  wake 
versus  downstream  distance.  Wake  width  at  a given  downstream  station  is 
defined  as  the  distance  from  the  axis  of  symmetry  at  which  the  cross- 
wake profile  of  axial  turbulent  intensity,  z » falls  to  half  of  its 

peak  value.  The  predicted  wake  width  at  z/R^  ■ 12  differs  by  only 
about  5%  from  the  value  presented  by  Lin,  et  al.  (1974). 

4.4  Calculation  Initialized  By  Boundary-Layer /Propeller  Model, 

Re  - 1.5  x IQ7 

The  Reynolds  number  for  the  calculation  described  in  this 

section  was  achieved  by  multiplying  all  lengths  in  the  calculation  for 
4 

Re  « 1.5  x 10  by  100,  and  all  velocities  by  10.  The  grid  was  again 

60  x 32,  but  since  high  resolution  near  the  body  was  required,  the  East 

boundary  was  set  to  z/Rb  * 12  and  the  North  boundary  to  z/R^  ■ 1.2  . 

There  were  20  points  between  z/Rj^  ■ 0.  and  0.6  , 24  points  between 

r/R^  • 0.  and  0.72  . Because  of  the  higher  space  resolution,  the  time 

step  required  for  convergence  was  one-fifth  of  that  used  in  the  calcula- 

4 

tion  for  Re  - 1.5  * 10  . All  other  aspects  of  the  calculation  for 

7 4 

Re  ■ 1.5  x 10  were  the  same  as  for  Re  ■ 1.5  x 10  . The  calculation 

required  about  55  minutes  of  CDC  7600  central  processor  time. 

A plot  of  the  calculated  wake  width  for  Re  ■ 1.5  x 10  ^ is 

4 

presented  in  Figure  4.12,  alongside  the  wake  width  plot  for  Re  ■ 1.5  x 10  . 

As  one  might  expect,  the  wake  is  narrower  for  the  high  Re  case. 


"P — 


Flow  Research  Report  No.  69 
December  1978 


5.  Conclusions 

A two-part  procedure  for  modeling  the  turbulent  flow  in  the 
near-wake  of  a submerged,  axlsymmetrlc,  self-propelled  vehicle  has  been 
presented.  Only  gross  parameters  such  as  the  body  speed,  the  body 
dimensions,  and  the  propeller  geometry  are  required  to  Initialize  the 
model.  The  first  part  of  the  procedure  establishes  conditions  at  the 
tall  of  the  vehicle,  just  aft  of  the  propeller,  which  can  be  used  to 
Initialize  the  second  part  of  the  procedure,  a finite  difference  simu- 
lation of  the  turbulent  near-wake.  Laboratory  data  can  also  be  used  to 
initialize  the  wake  simulation. 

A computer  code  based  on  the  procedure  has  been  used  to 
generate  four  sets  of  predictions,  three  of  which  have  been  compared  to 
laboratory  measurements.  TVo  sets  of  these  predictions  were  also 
Initialize!  by  the  laboratory  data.  As  can  be  seen  from  Sections  4.1, 
4.2,  and  4.3,  the  comparlsions  suggest  that  the  model  holds  substantial 
promise  of  being  able  to  successfully  describe  turbulent  near-wakes. 

In  spite  of  these  good  results,  several  aspects  of  the  model 
clearly  require  further  consideration.  For  example,  the  analysis  used 
to  construct  the  first  part  of  the  modeling  procedure  is  relatively 
simplified  and  is  only  justified  if  it  turns  out  that  the  flow  in  the 
turbulent  wake  does  not  depend  too  strongly  on  the  details  of  the  flow 
at  the  tail  of  the  vehicle.  The  degree  of  this  dependence  has  to  be 
determined  by  systematically  evaluating  the  sensitivity  of  the  down- 
stream flow  to  the  conditions  at  the  vehicle's  tail.  If  it  turns  out 
that  certain  upstream  variables  are  critical  for  the  wake's  evolution, 
then  these  quantities  must  be  more  accurately  determined.  It  is  our 
conjecture  that  the  mean  flow  quantities  are  far  more  critical  than  the 
turbulence  variables  for  the  wake's  evolution,  and  thus,  that  our  model 
would  benefit  from  a more  accurate  determination  of  the  mean  flow  just 
aft  of  the  propeller. 

A related  question  is  how  critical  are  the  methods  used  to 
complete  an  incomplete  set  of  laboratory  data  which  is  used  to  initialize 
the  calculation  of  the  turbulent  near-wake.  For  example,  in  the  simu- 
lation of  Schetz's  flow  we  have  used  an  eddy  viscosity  assumption  to 
determine  the  initial  R^  profile,  even  though  we  know  that  for  com- 
plex turbulent  flows,  eddy  viscosity  assumptions  often  fall.  Also,  in 
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Schetz's  case,  to  determine  the  Initial  dissipation  rate  profile,  we 
have  used  the  high  Reynold's  number  assumption  that  e ■ k /£  , where 
k Is  the  turbulent  kinetic  energy  and  £ Is  a turbulence  integral 
scale,  whereas  In  the  case  of  Chevray's  flow,  we  have  followed  Pope  and 
Whitelaw  (1976)  and  assumed  that  the  dissipation  rate  is  equal  to  the 
production  rate  of  turbulent  kinetic  energy.  It  is  not  clear  which 
approach  for  estimating  e is  better.  A final  example  is  that  we  have 
set  the  initial  mean  radial  velocity  profile  to  zero  for  Schetz's  flow, 
even  though  Pope  and  Whitelaw  (1976)  show  that  the  downstream  evolution 
of  Chevray's  flow  is  fairly  sensitive  to  this  quantity.  Our  hope  is 
that  the  downstream  character  of  the  wake  does  not  depend  strongly  on 
crudely  estimated  quantities  such  as  R^q,  E»  and  U in  the  above  examples. 
Sensitivity  studies  are  again  called  for. 

In  the  near-wake  calculation,  we  assumed  that  axial  gradients 
of  turbulence  quantities  which  appear  in  turbulence  diffusion  terms  may 
be  neglected.  Although  the  good  agreement  between  our  calculation  and 
Pope  and  Whitelaw' s for  the  Chevray  case,  suggests  that  this  approximation 
is  a good  one,  it  should  be  more  carefully  evaluated.  Further  numerical 
studies  should  provide  some  insight  into  the  applicability  of  and  the 
error  resulting  from  this  assumption. 

Finally,  we  come  to  the  seven-equation  turbulence  model  that 
is  used  in  the  near-wake  calculation.  Although  the  model  we  have  used 
here  is  derived  from  a large  number  of  flow  simulations,  it  would  be 
overly  optimistic  for  us  to  expect  that  it  is  necessarily  applicable  to 
the  flows  of  interest  here.  We  should  expect  that  the  model  would 
benefit  from  at  least  some  modification  or  "tuning"  of  constants. 

Because  of  the  short  time  span  of  this  effort,  none  was  attempted. 

Perhaps  after  applying  the  model  to  other  flow  situations  similar  to 
those  of  interest  here  (with  the  expectation  of  encountering  difficulty), 
we  will  acquire  the  necessary  insight  and  experience  required  to  make 
rational  and  helpful  model  improvements. 
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Figure  1-1.  Submerged  Self-Propelled  Body 


(3)  Repidly  Varying  Near  Wake. 

Figure  1-2.  Region  IK  • Key  Elements. 
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Fig.  2.1  Correlation  of  Aft-body  Boundary  Layer  Data  with 
Flat-Plate  Boundary  Layer  Approximations. 
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Aft-Body  Boundary  Layer. 
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Figure  4-6.  Tangential  Shear  Stress  Across  Self-Propelled  Body 
Wake  at  z/R^  « 10. 
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Figure  4-8a.  Radial  Shear  Street  Across  Prolate  Spheroid 
Wake  atz/Rb  - 6. 
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Figure  44b.  Radial  Shear  Strata  Across  Prolate  Spheroid 
Wake  at  z/Rb  - 12. 


( 





Figure  4-12.  Calculated  Self-Propelled  Body  Wake  Width  Versus  Downstream 
Distance  for  Re  = 1.5  x 10*  and  Re  = 1.5  x 107.  An 
Experimentally  Obtained  Value  for  Wake  Width  at  z/Rf,  = 12. 
for  Re  = 1.5  x 10*.  is  Also  Presented. 
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Appendix  A.l.  Estimate  of  Turbulent  Boundary  Layer  Entrainment 

Even  the  most  simple  turbulent  boundary  analysis  (see  H.  Schlichting' s 
Boundary-Layer  Theory,  McGraw-Hill,  1968,  p.  599)  results  in  an  expression 
for  the  boundary  layer  growth  d6/dz*  (where  z*  is  the  streamwise 
distance  from  the  leading  edge), 


with  C 


d6  „ 36  , 

dl*  T cf(z  } ’ 


» r-  (=  C'  in  Schlichting's  notation).  For  a non-flat 


surface  we  have  that  the  rate  of  change  of  the  boundary  layer's  volume 


flux,  \|>e  , is 


= wed6dy  _ d6 
dA  dz*dy  We  dz*  * 


where  dA  is  a differential  surface  area.  Thus 


Ct 

Vz*>  = f-  J we(5)  cf (O-^f  dC 


Now,  let's  assume  that  w (£)  ■ w (1  + ef(£))  £<<1  and  that 

0 00 

dA/d C “ KA(1  + e8(S))  (for  example,  on  a flat  plate  or  an  aligned 
cylinder  e = 0) , then 

z* 

V**)  " 7^  WeKA  / 


r weA(2*)  1 z*  / Cf<^ 


since  A(z*)  - KAZ**  Note  that  the  factor  in  the  brackets  is  the  definition 

of  the  total  friction  drag  acting  over  the  length  z*  (denoted  C^, 

here  and  in  Schlichting,  to  distinguish  it  from  our  or  Schlichting's 


\ 
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H>e  - 5.15  we(z*)A(z*)CF(z*) 

Because  of  the  crudeness  of  Schlichtlng's  analysis,  let's  discount  the 

accuracy  of  the  constant  5.15  but  still  retain  the  dependence  on 

w , A and  C_  . 
e F 

To  test  this  hypothesis,  we  will  use  Patel's  axisymmetric  6-to-l 

ellipsoid  data.  For  this  data,  < r*)  was  determined  from  integrating 

his  boundary-layer  velocity  profile  data  from  the  wall  location  to  where 

2_  2 

the  measured  stagnation  pressure  Cp  - [(P-Poo)  + 4pw  J /J spw^  is 
(0.995)^  (i.e.,  the  velocity  ratio  is  0.995  if  the  flow  were  on  a flat 
plate).  That  is, 

I'  w(y) 

J 2ir(r*  + y)  dy 

oo  woo 

0 

and  where  r*  is  the  body  radius  measured  normal  to  the  local  surface. 

To  estimate  Cp(z*)  utilize  Patel's  survey  station  amidship  (where 
the  boundary-layer  thickness  is  very  small  compared  to  the  local  wall 
radius)  and 


C 


F 


2 


e (z*) 

z* 


where  6 is  the  measured  momentum  thickness  of  the  boundary  layer. 


•-/ 


In  doing  this  for  the  survey  station  at  z*/L  , - 4.50  x 10  , 

which  according  to  Schlichtlng's  Figure  21.2,  which  shows  friction  drag 
measured  as  a function  of  Reynolds  number  based  on  flat-plate  length, 
gives  a value  of  wooz*/v  - 1.2  x 10^  which  is  50  percent  greater  than 
the  actual  value  of  w^zVv  . This  difference  is  presumably  due  to 
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Patel's  boundary  layer  trip.  To  compute  the  value  of  C at  other 

r 

stations,  use 


<*  *> 


effective 


1.2  x 10  (■ 


6 .1.71  + 


and  then  look  up  on  the  figure  what  (C  ) should  be  at  that  station. 

r r « r • 

With  the  preceeding  description,  the  following  table  of  values  was 
computed  for  Patel's  axisymmetric  boundary  layer. 


Survey 

Station 

Upstream 

Surface 

Area 

\p  /w 

Te  oo 

^F^F.P. 

ip  /w  AC 
ve  <*>  F 

(z*/L) 

(ft2) 

(ft2) 

— 

— 

0.662 

7.51 

0.116 

4.50  x 10-3 

3.43 

0.80 

9.11 

0.153 

4.35 

3.86 

0.85 

9.58 

0.158 

4.31 

3.83 

0.90 

9.98 

! 0.181 

4.28 

4.24 

0.93 

10.16 

0.204 

4.25 

4.72 

0.96 

10.32 

0.205 

4.23 

4.50 

0.99 

10.75 

0.206 

4.22 

4.55 

The  average  of  the  values  shown  in  the  most  right-hand  column  is  4.16 
with  a standard  deviation  of  11  percent  of  the  average.  To  see  how  this 
value,  4.16,  compares  to  flat-plate  turbulent  boundary- layer  data,  use 
Cole's  "universal"  turbulent  boundary  layer  profile 


. I ln  _I  + 5 + M 


where 


K - 0.41  (Karman's  constant) 


I 
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tt  = 0.55  (Cole's  pressure  grad,  parameter 
for  a flat-plate  B.L.) 

*■  T^/p  (The  "Friction  Velocity") 

1 2 

Note:  **  2(uT/we)  . The  quantity  ij^/w^LCj,  has  been  evaluated 

using  this  equation  in  a manner  quite  similar  to  that  used  for  Patel's 
data.  The  results  of  this  calculation  are 

1 LC 

Cf  ReL  0/6  = *e/veS  ^e/LCFWe 

3.2  x 10-3  1.5  x 106  0.106  0.793  3.73 

j 2 x 10~3  107  0.094  0.785  4.19 

Again,  a value  of  about  4 is  obtained  for  the  ratio  ip  /w  AC_ 

e e F 
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Appendix  A. 2.  The  Flow  Field  Induced  by  a Marine 

Propeller  Immersed  Within  a Rotational 
and  Nearly  Inviscid  Flow  Field 

In  this  appendix,  an  exact  linearized  solution  for  the  nearly 
inviscid  flow  through  a marine  propeller  is  derived.  This  solution  does 
rely  upon  certain  assumptions  in  order  to  make  the  problem  tractable. 
These  assumptions  are: 

1.  The  net  thrust  of  the  propeller  is  small  enough  so  that  the 
contraction  of  the  stream  tubes  passing  through  the  propeller 
can  be  neglected. 

2.  The  propeller  has  a large  number  of  blades  so  that  circum- 
ferential derivatives  can  be  deleted  from  the  equations  of 
motion  and  so  that  the  vortex  sheets  shed  downstream  from  the 
blades  can  be  approximated  by  a continuous  vorticity  distri- 
bution. 

3.  The  fluid  flow  is  inviscid  everywhere  except  at  the  plane  of 
the  propeller  where  viscosity  is  needed  in  order  to  invoke  the 
Kutta  condition  at  the  blade  trailing  edges  so  that  a net 
thrust  and  torque  act  on  the  fluid.  Also,  viscosity  has  acted 
upstream  of  the  propeller  in  order  to  create  the  rotational 
flow  approaching  the  propeller. 

4.  The  flow  field  is  axially  symmetric  with  no  forebody  present 
near  the  propeller  plane. 

A cylindrical  polar  coordinate  system  (z,r,0)  is  chosen  with  its 
origin  at  the  intersection  of  the  axis  and  the  propeller  plane  and  with 
positive  z measured  downstream  of  the  propeller;  w,u,v  are  the  fluid 
velocity  components  in  the  three  coordinate  directions;  £,n,5  are  the 
components  of  the  fluid  vorticity  vector  along  the  coordinate  directions. 
The  Inviscid,  axially-symmetric  equations  of  motion  for  an  incompressible 
fluid  are  then; 


(Continuity) : 


h (wr>  + h (ur)  ‘ 0 
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(Axial  Momentum) : 

(Radial  Momentum) : 

(Angular  Momentum): 
(Streamwise  Vorticity): 
(Radial  Vorticity): 
(Circumferential  Vorticity): 


dz  dr  r 
3v  , 3v  , uv  . 

“ 9?  + u 3?  + r • 0 
5 i|^  + - , 

dr  r 

_ - 3v 
n " 3z  ’ 


3u  3w 
3z  3r 


In  the  above,  all  length  scales  (z  and  r)  are  made  dimensionless  by 
dividing  by  the  propeller  radius  R and  all  velocities  are  made  dimen- 
sionless by  dividing  by  the  freestream  speed  wm  . The  continuity 
equation  is  satisfied  by  introducing  the  perturbation  stream  function 
ip(z,r)  defined  by: 


w(z,r) 


wo(r) 


1 3£ 

r 3r 


u(z,r) 


1 3£ 

r 3z 


(A-l) 


where  wQ(r)  is  the  axial  velocity  profile  far  upstream  of  the  propeller. 
Substituting  Equation  (Al)  into  the  circumferential  vorticity  equation 
results  in 


dz2  V 3r  'r  3r; 


(A-2) 


Cross-differentiating  the  axial  and  radial  momentum  equations  in  order 
to  eliminate  the  pressure  field  results  in 


w + u - (|£-  + uc)  - 0 

dz  3r  r '3z  w » 


(A-3) 
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whereas  the  angular  momentum  equation  can  be  rewritten  as 


w K + u 1£  + i ill  . _n  3w  ; 

W 3z  U 3r  r 3r  0 3r  0 


(A-4) 


9w 

(n^—  vanishes  upstream  and  is  near  zero  downstream  of  the 
dr 

If  the  propeller  thrust  is  small  we  expect  that  everywhere 
w ~ 0(1)  u,v  <<  0(1)  . Thus,  we  can  approximate  £ and 


Equations  (A-3)  and  (A-4)  as 


propeller) . 
in  the  fluid 
£ using 


«o(r> 

So(r)  + Cj_(r) 


£(z,r)  = 


r 

(q(r) 


z < 0 
z > 0 


z < 0 
z > 0 


(A- 5) 


where  is  the  initial  circumferential  vorticity  distribution 

dw  ' 

(£q  ■ - ) and  C^(r)  and  £^(r)  are  the  changes  caused  to  the 

circumferential  and  streamwise  vorticity  (the  latter  is  related  to  the 
swirl  velocity)  by  the  propeller.  Since  the  present  solution  is  only 
intended  to  be  applicable  at  and  upstream  of  the  propeller  plane,  the 
parallel  flow  approximation  implicit  in  Equation  (A-4)  is  not  as  restrictive 
as  it  might  first  seem.  Since  the  upstream  flow  only  changes  because  of 
the  downstream  vorticity  field  (i.e.,  Biot-Savart's  law  of  the  induced 
velocity  from  a vortex  line)  Equation  (A-4)  is  tantamount  to  assuming 
the  induced  velocities  upstream  are  nearly  the  same  as  they  would  be  for 
a parallel  flow.  Strictly  speaking,  a parallel  flow  would  require  that 
the  radial  velocity  be  zero  everywhere,  whereas  a finite  radial  velocity 
distribution  is  still  calculable  from  the  present  analysis.  If  Equation 
(A-5)  and  (A-2)  are  now  combined  we  arrive  at 

72  $ + h <¥  U>  ■ -hM  H">  • <A-6) 

3 z 

where 


( 
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H(z) 


0 

1 


2 < 0 
2 > 0 


The  perturbation  streamfunction  i|j(z,r)  can  now  be  separated  into  the 
sum  of  two  functions,  one  of  which  is  independent  of  z and  the  other 
of  which  is  an  odd  function  of  z . Far  upstream  this  sum  vanishes 
(i.e.,  there  is  no  perturbation  far  upstream)  and  far  downstream  the  two 
add  up  to  represent  the  jet-like  flow  caused  by  the  propeller.  The 
axially  independent  portion  of  this  perturbation  streamfunction  is 
related  to  the  induced  axial  velocity  at  the  propeller  plane  (which  is 
one-half  the  axial  velocity  at  z = +»  if  the  flow  were  always  inviscid) 
through 


^(r) 


00 

- if 


(p)pdp 


(A-7) 


where  Aw+(p)  is  the  far  downstream  axial  velocity  distribution  for 
a completely  inviscid  flow. 

The  odd  contribution  to  the  perturbation  streamfunction  may  be 
shown  to  be  given  by 


00 

f 


+Xz . 


</»2(z,r)  - + ^(r)  + r I f(X)e  J^XrJdX  , z > 0 


(A-8) 


with 


f(X) 


oo 

- if  4»+ 


(p)JQ(Xp)pdp  , 


from  a separation-of-variables  solution  to  the  resulting  partial  dif- 
ferential equation;  Jq  and  are  Bessel  functions  of  the 

first  kind  of  order  0 and  1. 

The  induced  velocity  field  upstream  of  -he  propeller  can  then  be 
determined  from  the  perturbation  streamfunction  if  the  (hypothetical) 
axial  velocity  profile  Aw+  can  be  specified.  The  upstream  axial  and 
radial  velocities  perturbations  are: 


( 
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00  CO 

Aw(z,r)  m\J  Aw+(p)pdp  J XeXzJQ(Xr)Jo(Xp)dX  , 
o o 

oo  oo 

u(z,r)  = - Aw+(p)pdp  j XeXzJQ(Xr)J1(Xp)dX  . 
o o 

With  a considerable  amount  of  work,  these  relations  can  be  shown  to  be 
equal  to 

OO 


(A-9) 


Aw(z,r ) m f Aw+(p)  Gw(z,r;p)pdp 
o 

oo 

- 2i rf  4"+ 


z < 0 


(A-10) 


u(z, r) 


(p)  Gu(z,r;p)pdp 


where 


Gw(z,r;p)  *=  E^m>  , ~~ 

z2  + (r  - p)2\Jl2  + (r  + '2 


P)‘ 


Gu(z,r;p)  = (r2  * p2  - z2)Gu(z,r;p)  + 


K(m) 


JJl  (r  + p)2 


and  where  K(m)  and  E(m)  are  complete  elliptic  integrals  of  the  first 

2 2 

and  second  kinds  with  modulus  m = 4rp/(z  + (r  + p)  ) . (See,  for 
example,  "Handbook  of  Mathematical  Functions"  by  Abramowitz  and  Stegun.) 
At  the  propeller  plane,  z = 0 , these  simplify  somewhat  to 


Aw(o,r)  « j Aw+(r)  , 

OO 

u(°.r>--fe/  Aw(0*p)  [f^p  + 7^]  pdp  • 


(A-ll) 


m ■ 4rp/(r  + p)z  . 


In  other  words,  the  Induced  axial  velocity  at  the  plane  of  the  propeller 
is  one-half  the  (hypothetical)  far  downstream  axial  velocity  and  the 
induced  radial  velocity  is  directly  related  to  the  induced  axial  velocity 
at  the  propeller  plane.  Note  also  that  nothing  has  as  yet  been  said 
about  the  swirl  velocity  distribution  v(z,r)  . From  Equation  (A-5) 
along  with  the  definition  of  the  streamwise  vorticity  it  follows 


^ 5TT- 
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v(z,r)  = v+(r)  H(z)  . 


(A-12) 


In  other  words,  once  generated,  the  swirl  velocity  profile  is  independent 
of  z (again,  subject  to  the  approximation  that  the  flow  is  always 
inviscid) . There  are  no  net  induced  velocities  upstream  of  the  propeller 
caused  by  the  swirl  velocity  vortex  lines. 

Note,  however,  that  the  swirl  velocity  distribution  is  related 
to  the  axial  velocity  perturbation  through  the  parameters  which  characterize 
the  propeller  and  its  blades.  In  fact,  if  these  characteristics  are 
specified,  it  is  possible  to  utilize  the  present  formulae  to  compute  the 
axial,  radial  and  swirl  velocity  fields  as  well  as  the  total  thrust  and 
torque  acting  on  the  fluid  by  a procedure  quite  similar  to  that  described 
by  Schwartz  and  Bernstein  (1975).  They  addressed  the  propeller  flow 
field  utilizing  two-dimensional  airfoil  theory,  and  developed  a computer 
code  based  on  their  analysis.  As  mentioned  in  the  text,  the  above 
analysis  is  in  some  ways  more  complete  than  the  previous  work,  and  a 
modification  or  a rewrite  of  the  Schwartz-Bernstein  code  was  considered. 

It  was  concluded,  however,  that  such  an  effort  would  result  in  only  a 
slight  improvement  over  the  axial  and  swirl  velocities  currently  computed. 

On  the  other  hand,  the  existing  code  dotB  not  compute  the  radial  velocity 
field  at  the  plane  of  the  propeller,  but  assumes  it  to  be  identically 
zero.  A better  estimate  of  this  prof ile  is  available  from  Equation  (A- 
11). 

To  illustrate  tl.e  magnitude  and  shape  of  the  radial  velocity 
distribution  an  exact  evaluation  of  Equation  (A-ll)  is  possible  for 


certain  choices  of  Aw(o,r)  . For  example,  if  Aw(o,r) 
radial  velocity  distribution  can  be  shown  to  be: 


'1  - r 


u(o.r)  _ 

Aw 

max 


-irr/8  , 0 < r < 1 


— F(—  — • — • — ) r > 1 

2 v2*  2*  29  2'  9 L d- 
6r  r 


( 
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where  F is  a hypergeometric  function.  Thus,  the  radial  velocity  is 

inward  and  varies  linearly  over  the  face  of  the  propeller.  At  a large 

_2 

distance  from  the  axis  the  radial  velocity  approaches  zero  like  r 
It  is  noteworthy  that  the  maximum  radial  inflow  velocity,  which  occurs 
for  r slightly  greater  than  1,  is  slightly  in  excess  of  ir/8  times 
the  change  in  the  axial  velocity  which  illustrates  that  the  changes  in 
either  the  radial  or  axial  velocity  are  of  the  same  order.  A numerical 
calculation  of  the  upstream  induced  velocities  (from  which  the  pressure 
field  can  be  computed)  has  been  obtained  which  verifies  that  the  upstream 
influence  is  only  appreciable  for  one  propeller  diameter  (i.e., 

|C  I < 0.01  for  z < -2). 

'p' 
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Appendix  A. 3.  Derivation  of  Jump  Equations  Across  Propeller  Plane 
for  Grabowski's  Turbulence  Equations 


Assume  that  the  axial  flow  speed,  W(z,r)  , and  the  radial  flow 
speed,  U(z,r)  , are  continuous  functions  across  the  propeller  plane 
located  at  z=0  whereas  the  swirl  flow  component  V(z,r)  goes  from 
zero  ahead  of  the  propeller  to  V(r)  at  z = 0+  . Now,  assume  that 
each  of  the  turbulent  velocity  correlations  on  any  streamline  R 


R 


zz 


R , Rq 
rr  0z 


ee  » 

Rgr  , and  Rzr  can  experience  a step  change  as  the 


streamline  passes  through  or  near  the  propeller  plane  and  similarly  for 
the  dissipation  e . 

Consider  now  the  P.D.E.'s  for  the  seven  turbulent  quantities 

mentioned  above  which  are  shown  on  page  76  of  the  ICWAKE  report.  Note 

the  left-hand-side  (LHS)  of  these  equations  all  involve  the  material 

9 9 

derivative  operator  U -r  W . Using  the  continuity  equation, 

this  operator  can  be  written  (using  the  equation  for  R^  as  an  example) 
as 

~ .UR 

4—  (UR  ) + (WR  ) — . 

9r  rr'  9zv  rr  r 

If  this  operator  is  integrated  from  z = -A  to  z = +A  we  have 

A A 

4-  / UR  dz  + [wR  1 - - / UR  dz 

dr  J rr  I rr  J r J rr 

-A  -A 

denotes  the  jump  in  the  enclosed  quantity  across  the  propeller 


where 

plane.”  This,  in  turn,  can  be  written  as 

/ 

A 

d I 1 
r 

ar  ) r 

-A 


/“r r4*!  + [“\r]  ; 


— U (0,r)  (R-  + R )A  + - U(0,r)R-  A 
r rr  rr  “ r * rru 


•r)  [Rrr]  + r h 

;»(0.r)[Rtr]  + |i  U(0,t)  [«„]  ) A 
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Thus,  by  allowing  A to  approach  zero,  the  second  term  In  the  last 

relation  will  vanish  whereas  the  first  term  remains  finite.  In  this 

manner,  the  LHS  of  each  of  the  seven  P.D.E.'s  when  integrated  across 

the  propeller  plane  reduces  to  W(0,r)  times  the  jump  in  the  quantity 

3 8 

being  operated  on  by  U + W . 

Now,  consider  the  right-hand-sides  of  the  seven  P.D.E.'s  when 
integrated  across  the  propeller.  The  analysis  for  a typical  production 
term  is 


A 

/ Prr(z,r)dz 
-A 


-2 


/ (R, 


3U  3U 

3r  rz  3z 


-2(R~  + f R 1 )4p 

rr  rr  3r 


z=+A 


A + 2R~ 

rr  3r 


z=-A 


2 (R  + 
rz 


["r  1 

L rzJ 

’ 3z 

z«+A 


4 + 2R-  |2 

rr  dz 


z=-A 


+ 2(r;6 


+ 


A 


- 0(A)  . 

Here,  even  if  3U/3r  is  not  continuous,  all  terms  are  still  multiplied 
by  A which  can  be  made  as  small  as  necessary.  On  the  other  hand,  the 
analysis  for  another  typical  production  term  is 

A A 

/ W*  • -2  / 

-A  -A 


R 3V  3V  U 

r0  3r  R0z  3z  R00  r 


dz 


( 
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c 

; -2<N  + [ 

r 1 

r0j  dr 

- 4 - 2<R9; 

_ 2 U(0,r,) 
r 

J (R00  + 

1 

<] 

<D 

<D 

l i 

I 

■ - 2 <N + 

M > 

[v]  + 0(A) 

• *2  MM 

since  Rn  «= 

0z 

66 


0 for  z < 0 . 


Results  for  other  production  terms,  after  analysis  similar  to  the  above 
is  as  follows: 

A 

I 


P dz  = 0 
zz 


-A 


I Pr0dz  = -(Rrz  + 


-A 


M]  ,H 


f P dz  * 0 
J rz 


-A 


/ p0zdz  - “(Rzz 


-A 


+ M [v] 


/ 


P£dz  - -C. 


el^  CR0zMR0j  ) [v]  ^C£l^3 


M 


Ld 


where  k2  - \ (R^  + R00  + R^)  . 


(since  R~  » 0) 

wZ 


+ jm 


[Rez]  M 
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For  the  first  part  of  the  modeled  pressure-strain  correlations,  we  have 

A _ 

/ ^ij , ldz  ; 0 to  order  A , 

-A 

and  for  the  second  part  of  the  modeled  pressure-strain  correlations,  the 
following  results  hold: 


i w'-  /iv-j*  nw 

-A 


<sl"“  Pk  ’ T<Prr  + PO0  + P«>  and  / Drrd*  ‘ °«» 


c3  + c5 


M M 


r ' 2C3 

f 2C 

L. 

r 

1 ^QQ,2dz  ' - — 
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i 

+ 0(A) 
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■t'NM  -tW[v] 

4C  - 2C5  [- 

/ /¥*♦* ,[■*][»]  /¥- 


A ^ 2C  A P 

f W ^ /- 


2Cc  “ Pc 


* 3 <2C5  " C3> 


Kl[v]  +2s  W[v]  - £[■*][*] 

- V [»ej  ] 


/ *re,2d*  ■ C3(llrz  + [Rrz]  >[V] 
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/ *rz,2iz  ' C5  M[V 


/ W"2  ’ C3(RIz  + [Rzz]  > [V]  - C4((k'>2  +[1‘2])  [V] 


+ c5(R00  + 


<C3Rzz  ' C4<k">2  + S 


+ [vMv] 

W [v]  ■ 


Using  the  above  results,  ignoring  j^R^J  , R0eJ  , RZz j ’ [ £j  ’ j^Rrz' 
relative  to  R , Rfifl  , R , e , R , and  ignoring  the  diffusive 

2T1T  v/ v ZZ  Tl  Z 

A 

terms  (i.e.,  the  J Df rr.dz  etc.)  , we  have  the  following  new  forms  for 


the  seven  turbulence  equations  after  they  are  integrated  across  the 
propeller  plane. 


* M 

“ [Ree] 

“ M 
“ W 

“ m 


!<C3  + C5>  [R8z]  [V] 

2 M [']  + ^ 


4C3  - 2S  r in  4S  - 2C5  - 6 

■V  WH  - 


MM 


\ «c5  - C3>  [Rez]  [V] 


-R  V + C,R 
rz  I 3 rz 
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[v]  - <c3  - i>r;z  [v] 
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M [«»z]  ■ -RL  [ v]  + <c3*;,  - Vk'>2  + SRee>  [v] 


+ <c5  - r)Ree 

L.  -J 

-^R- 
2 rr 

V 

1 

F(r)  V , where  r « (r  Is  the  propeller  radius), 
L J P P 

and  lastly, 

“ H ' - [v  ] 


From  the  previous  equations, and  the  constants  given  in  Appendix  B (p.  90), 
the  following  conditions,  presented  in  Section  2 (p.  16),  may  be  obtained. 
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M 
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Appendix  B:  Turbulence  Model  Equations  For  Axisymmetric, 
Incompressible  Flow  With  Swirl 


The  turbulence  model  equations  presented  in  this  Appendix  are  based 
on  the  model  of  Hanjalic  and  Launder  (1972)  and  its  extension  by  Launder, 
Reece  and  Rodi  (1975) . The  equations  are  presented  in  the  general  form 

Advection  - Production  + Dissipation  + Pressure-Strain  (where  applicable) 

+ Diffusion 

Dissipation  Rate  Equation 


ul^+W^'P  - C.~+Df 
3r  3z  e e2  k £ 


Turbulent  Correlation  Equations 


U if  + “ IT  - 2 7 \e  ■ Prr  - 3 e + *rr,l  + »rr.2  + D£rr 

V 2 ~ ~ 

u “37“  + w 1 T + 2 7 Rr0  ‘ pee  " 3 e + *00,1  + *00,2  + Dfee 

3R  3R  0 

U + W - P - 4 e + 4)  ,+4>  _ + Df 

3r  3z  zz  3 zz,l  zz,2  zz 

3F  o 3R  A v / \ ~ ~ 

u ~TT  + w 17"  + 7 (Rrr  " Ree)  = pre  + *re,i  + *r0,2  + Dfre 


OR  OR  ..  « 

U + w - - R.  - P + <t>  . + <f>  » + Df 

3r  3z  r 0z  rz  rz,l  rz,2  rz 

3Ra  3Rfl  v „ 

U + W -r22-  + - R **  P„  + <t>ft  , + 4>0  0 + Dffi 

3r  3z  r rz  0z  0z,l  0z,2  0z 


The  production  terms  in  the  above  equations  are  defined  as 


P£  * - C£1 
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3z 
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3r  rz 

3r  6z  ( 3z  r / J . 

The  first  part  of  the  modeled  pressure-strain  correlations  was  suggested 
by  Rotta  (1951): 


+rr,l--Clt(Rrr-!k)  > 

*ee,i  "'if  (Ree  " I k)  • 
♦„,i  ■ - ci  t (\*  - 1 k)  • 
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<t>  ! 

rz,  1 


" C1  k Rrz  * 


and 


*6z,l  " " C1  k R6z  ’ 


The  second  part  of  the  modeled  pressure-strain  correlations  is  due  to  Launder, 
Reece  and  Rodi  (1975) : 
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K00 ,2 
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where  is  the  production  rate  of  turbulent  kinetic  energy 

Pk-!(PrI  + P8e  + Pzz)  • 


and  the  terms  are  defined  as  follows: 
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The  expression  for  the  diffusion  of  dissipation  rate  is  based  on  the  proposal 
of  Hanjalic  and  Launder  (1972)  : 

m=-ts [ i(" If + *„ H)] ♦ ce £ [!(»„ If ♦ It )]  • 


Flow  Research  Report  No.  69 
December  1978 


-87- 


When  n ■ 0 the  turbulent  correlation  diffusion  terms  below  are  based 
on  the  proposal  of  Daly  and  Harlow  (1970);  when  n * 1 , they  are  based  on 
the  proposal  of  Hanjalic  and  Launder  (1972). 
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There  are  six  model  constants  which  must  be  specified. 
Following  Hanjalic  and  Launder,  we  let 


"7 -r'r^ 
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C£  = 0.15  , 


C£1  - 1.44  , 


Ce2  = 1.90  , 


and  following  Launder,  Reece  and  Rodi,  we  let 


Cx  = 1.5  , 


C2  = 0.4  , 


C _ = 0.25  when  n = 0 
sO 


C , = 0.11  when  n = 1 
si 


The  constants  and  are  defined  as 


C3  = (C2  + 8) /II  , 


C4  = (30C2  - 2)/ 55  , 


C5  = (8C2  - 2)/ll  . 


In  order  to  simplify  this  rather  lengthy  system  of  equations,  we 
assume  that  in  the  turbulent  diffusion  terms,  axial  derivatives  of  the 
various  turbulent  quantities  may  be  neglected  compared  to  their  radial 
derivatives.  We  emphasize  that  this  assumption  is  applied  only  to  the 
turbulent  diffusion  terms. 

The  simplified  turbulent  diffusion  terms  are  as  follows: 
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Appendix  C:  Coordinate  Transformation  Parameters 

This  Appendix  describes  the  calculation  of  the  four  parameters 
appearing  in  the  logarithmic  coordinate  transformations  used  in  this 
study.  The  transformations  are  briefly  described  in  Section  2.4  of 
this  report. 

We  introduce  the  new  independent  variables  y and  x in  place 
of  r and  z where 

y = ay-1ln(l  + r/by) 

and  x = a 1ln(l  + (z  - Z )/b  ) . 

X 1 X 

In  order  to  arrive  at  appropriate  values  of  the  parameters  ay 

and  by  which  appear  in  the  r to  y transformation,  we  first  specify 

Rmax  location  of  the  North  boundary),  and  the  value  desired  for 

y at  R .We  will  denote  this  value  as  y .We  assume  that  a 
max  ■'max 

total  of  Mc  grid  points  will  be  distributed  radially  over  the  region 

0 < r < R where  R is  less  than  R . Generally,  R is  approxi- 

mately  the  wake  radius  at  the  West  (upstream)  boundary,  and  Mc  is  the 

number  of  grid  points  across  the  wake  at  this  station.  Then,  since 

Ay  * y / (M - 1)  , y , which  is  the  y value  corresponding  to  R , 
max  c c c 

is  given  by  yc  * (Mc  - l)Ay  . After  inverting  the  transformation  we 
may,  therefore,  write 


R ■ b (exp(a  y ) - 1) 

c y r yJ  c ' 

and  R»ax  ' by(exp(Vn»ax)  “ * 

We  rearrange  and  divide  these  expressions  to  yield 
(1  - expta^))  Rc 
« - *"%“ ' ' 0 ' 

We  solve  this  equation  for  ay  using  a simple  Newton-Raphson 
iterative  procedure  beginning  with  an  initial  estimate  (guess)  for  a . 
With  ay  known,  we  obtain  by  explicitly  from  either 
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b “ R (exp(a  y ) - 1)  , 

y c r yJc  ’ * 

or  b ■ R (exp (a  y ) - 1) 

y max  r yJmax 

The  calculation  of  the  parameters  and  bx  of  the  z 

transformation  follows  a procedure  Identical  to  that  which  has 
described  above. 
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Appendix  D:  Finite-Difference  Equations 


The  finite-difference  formulations  of  the  transformed  version  of 
the  mean-flow,  pressure,  turbulence  model,  and  East  boundary  condition 
equations  are  presented  in  this  Appendix.  The  finite-difference  grid 
system  is  shown  in  Figure  3.2  and  has  been  described  in  Section  3 of  this 
report. 

D. 1 Notation 

In  the  difference  equations  presented  below,  centered  second-order 
differences  are  used  to  approximate  most  spatial  derivatives,  and  are 
denoted  as  follows  (where  $ can  represent  U , V , W , Rrr  , etc.): 


and 


3x  ~ Doi^lj  = ^i+1, j " ^i-l,3)/2Ax  » 

3y  = DoJ*iJ  = ^i, j+1  " *i,J-l),2Ay  ’ 

Z D44D-i*iJ  7 (<,*i+l,j  + *1-1, J " 
d2<t>  „ 


3y 


2 8 D+jD-J*ij  = ^*i, j+1  + *i,j-l  " 2<^ij)/Ay  • 


while,  first-order  backward  differencing,  which  is  used  at  i ■ N to 
approximate  certain  derivatives  with  respect  to  x , is  denoted  as 


3£  ~ 

3x  ~ 


D <t> 
-i9ij 


The  weighted  centered-upwind  finite-difference  formulation  used  to 
approximate  axial  convection  terms,  as  described  in  the  body  of  this 
report,  is  written  using  the  above  notation  as 


W 3x  * (WijDoi  ” “JWijl  2*  D+iD-i)*iJ  • 

This  formulation  makes  clear  the  introduction  of  artificial  diffusion 
for  non-zero  values  of  ; the  artificial  viscosity  at  grid  point 
i.J  is  aJw^lAx/2  . 

As  a result  of  their  definition  on  a staggered  grid  system  (with 
grid  points  denoted  as  i+*j,J-Mj  ) the  x and  y derivatives 
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of  P at  grid  point  i,j  are  approximated  by  averaged  centered  dif- 
ferences so  that 


dx  = (DoiP)ij  = [<Pi-Hs,  j-Hj  " Pi-l|tJ-U|> 

+ ^i+is.j-^  ' /2AX  ’ 

and 

If  * (DojP)ij  E [(Pi4*s,j4*s  " Pi*»*l,j->l) 

+ (Pi-»S,j-^  " Pi-»S,3-JS)]  /2Ay  ‘ 


Finally,  the  values  of  the  coefficients  f(y)  , g(y)  , h(y)  , 
s (x)  , and  t(x)  are  denoted  at  the  grid  points  by  fj  , g^  , h^  , 

8^  , and  t^  , respectively. 

D.2  Mean-Flow  Equations 

As  discussed  in  Section  3.2.2  of  this  report,  the  ADI  procedure  intro- 
duces two  half  time-steps  each  of  length  1/2  At.  During  alternate  steps, 
appropriate  y and  x derivatives  of  mean-flow  quantities  in  Equations 
(3.13),  (3.14)  and  (3.15)  are  differenced  implicitly.  The  implicit 
differencing  results  in  sequences  of  tridlagonal  matrix  equations  for 
Uij  * Vij  * an<*  Wij  which  are  easily  solved.  Note  that  in  the 
equations  presented  below,  the  Reynolds  stress  terms  are  approximated 
at  the  appropriate  time  level  and  at  grid  point  i,j  by 


T0M”  ' •iDol*rl  + Fj  Vt  + V'rr  * V • 


T0MV  ' 'AA.  + V»  + V*  ■ 


’ 'AA.  + fJDoA,  + hJBr*  • 


In  the  above  difference  expression  and  those  which  follow,  the  i,j 
subscripts  ere  omitted. 
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D. 2.1  ADI  First  Half  Time-Step 

During  the  first  half  time-step,  which  advances  the  solution  from 
time  level  n to  n+*s  , y derivatives  of  the  mean  velocities  U , 

V , and  W are  treated  implicitly.  The  following  finite-difference 
equations  are  thus  obtained: 


""wi "■ + v'V*'1 + D+iD-i>u"  - hj<v”)2 

- -VV>" + Re'1(  [ fjV-j  + + W-i  - h32J 

+ (s^D+iD_i  + t±Doi)Un  J - TURBU"  , 

V"m/2  V"  + £jra>ojV”1*!  + D«D-i)V"  + y°V” 

■ »*‘1  | [ fiV-i  + (8J  + fJhJ)DoJ  - hi]  ^ 

+ (sJd^D^  + t1Dol)Vn  - TURBVn  , 

and 

* fju°  + - 0il“”lr  D«D-1>H° 

- -i(0olp)”  + Re'1  ( [ fjV-J  + <*J  + fJhJ)DoJ  ** 

+ (sV.D  . + t.D  .)w"  - TURBW"  . 

1 +1  —1  1 Ol 

As  a result  of  the  use  of  implicit  differencing  in  the  y direction 
only,  these  three  difference  equations  may  be  rearranged  as  independent 
systems  to  the  form 

Aij*i,J-l  + + Cij*i,J+l  " DiJ  * 

where  2 < J < M - 1 , for  each  value  of  1,  2 < 1 < N - 1 . The 

quantity  ^ represents  or  , and  the 
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coefficients  and  D^  are  evaluated  using  the  known 

values  of  the  solution  at  time  level  n . It  should  be  noted  that  , 

Bij  * Cij  and  Dij  differ  depending  on  whether  represents  , 

or  . The  boundary  conditions  at  j+1  and  j - M are 

incorporated,  and  the  equation  systems  are  solved  suing  the  Thomas 
algorithm  (see  Roache,  1972)  at  each  value  of  i for  <J>^  where  1 < J < M . 

D.2.2  ADI  Second  Half  Time-Step 

During  the  second  half  time-step  which  advances  the  solution  from 
time  level  n-f^s  to  n+1  , x derivatives  of  U , V , and  W 
are  treated  implicitly.  The  following  equations  result: 


,n+l  „n+>s 


x un+1  - h (V®4**)2 

- + Ite'1  [ f jD+jD-j  + + WDoj  - 

x Un+Ss  + (s2D+iD_1  + tiDol)Un+1  - TURBUn44s  , 


vn+1  - vn+Js 


ATT r-  + + 8i^Doi  - V>-i> 


vn+1  + hjUn4V** 


• t*'1  + <*i + w\i  - h!  ]'?* 

+ (bJd+1D_1  + t1Dol)Vn  - TURBVn4i*  , 


W**1  - if* 
At7T 
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w"  + 1 
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As  is  the  case  during  the  first  half-step,  these  equations  can 
also  be  rearranged  as  Independent  tridiagonal  systems  but  of  the  form 


Aij*i-l,J  + Bij4>ij  + Cij*i+l,J 


i 5 ’ 


where  2 £ 1 < N - 1 , for  each  value  of  j,  2 £ j < M - 1 . Here 

, and  the  coefficients 


, . „n+l 

<p^  represents 


,,n+l 
• vij  or 


*ij  • Bij  • 


and  are  evaluated  using  the  known  values  of 


the  solution  at  time  level  n+*s  . When  the  boundary  conditions  at 
i - 1 and  i « N are  incorporated,  the  systems  are  solved  at  each 
value  of  j for  <f>^  where  2 _<  i < N (remember  that  and 

Wj  are  specified). 


D. 3 Turbulence  Model  Equations 

The  turbulence  model  equations  which  determine  the  seven  turbulence 
model  quantities  contain  second  derivatives  in  the  y direction.  They 
are,  therefore,  treated  implicitly  only  in  the  y direction,  and,  as  a 
result.  Identical  finite-difference  formulations  are  used  to  advance  the 
solution  from  time  step  n to  n-Mj  and  from  time  level  n-Hs  to 
n+1  . 


Since  the  difference  equations  for  the  seven  quantities  are  all 

quite  similar,  only  the  equation  for  R is  presented  here;  the  equa- 

rri1 

tion  is  written  as  follows 


e*.  r" 

rr  rr 

I72t 


n 

rr 


where 


rr 


- 

■PJr-f£°  + ^r,l  + *?.,2+DC  • 

- -2|ff.Rn  D + s.R  D .1  0"  - h.IBaVBi 

| L J rr  oj  i r*  oi J J r0  / 

n ' 


w ■ -ci  ^ (C  - ! k°>  • 
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^rr,2  ' ~C3^*rr  ' ! ^ “ Wo/  ~ S^r  “ ! *k>  * 


and 


D - -2f 
rr  j 


Rn  D ,U°  + RnftD  .Vn  + Rn  D .w” 
rr  oj  r8  oj  rz  oj 


The  finite-difference  approximation  to  the  turbulent  diffusion 


expression  Dfrr  , is 


Df  = (2n  + 1)C  h f D . 
rr  sn  j j oj 


h-1^ 
J c” 


f ,R  D ,Rn+4j  - 2h,(Rnft)2 
j rr  oj  rr  j'  r6 


- 2C  h — 
sn  j n 
J e 


(n  + 1) 


fvRnaD  .Rna  + h.R"  (Rn  - R* 
j r0  oj  r0  j 00*  rr  00 


+ n 


f.Rn  D .r"  + 2h.(R"  )2 

j rr  oj  ee  jv  r0' 


(The  quantity  n which  appears  as  a coefficient  in  the  above  expression 

does  not  refer  to  time  level.  It  is,  instead,  simply  a parameter  which 

may  take  on  the  values  zero  or  one.  See  the  discussion  of  the  turbulent 

diffusion  model  formulation  in  Appendix  B for  details.) 

At  i * N we  replace  the  weighted  centered-upwind  convection  by 

the  simple  first-order  expression  8..W...D  ,R  so  that  the  difference 

N Nj  -l  rrNj 

equations  are  defined  for  2 £ i £ N , 2 £ j £ M . 

The  difference  equation  for  Rrr  presented  above  (with  the  slight  modifi- 
cation at  i • N ) is  rearranged  into  independent  tridiagonal 

systems  of  the  form 


A..R  + B.,R  + C..R  - D..  , 

ij  rrltrl  ij  rr±J  ij  rri>j+1  ij 

where  2£J£M-1,  for  each  value  of  i , 2 £ i £ N . As  in  the 

case  for  the  mean-momentum  equations  during  the  first  half  time-step, 
these  equations  are  solved  at  each  1 for  1 £ j £ M after  the 
boundary  conditions  at  J ■ 1 and  j ■ M are  incorporated. 

The  difference  equations  for  the  other  six  turbulence  models 
quantities  are  similar  and  arc  solved  in  identical  fashion. 
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D.4  Poisson  Pressure  Equation 


The  transformed  Poisson  pressure  equation  (Equation  3. 16)  is  solved 
numerically  for  P^^  at  each  half  time-step  following  the 

calculation  of  both  the  mean  velocities  and  the  turbulence  quantities. 
The  finite-difference  formulation  of  the  pressure  equation  at  time  step 
n is  written  as  follows 


f j+*SD+jD-j  + (8j-^  + f j-HshWDoj  + 8i4isD+iD-i 
+ ti+isDoi]  Fi+*S,j+*S  Fi+J5,j+l5  ’ 


where 


Fi-H*,j+Js  “ fi4i5(DojCy)l+i5.j+J5  + (8j-^  + fj-^ 

X hJ+JS)Cyi4^f j44j  + 8i+*S(DoiCx)i4is,j4^ 

+ ^'l*.**"  ***'**  AT/2  ' 

* * 

The  Dqj  and  operators  have  been  defined  above;  and 


c n - -ke  n 

yi-ns,j44s  4 y±. 


“ +cn  +cn  +c“), 

i+l,j+l  yi+l,j  yi,j+l  yij 


Xi+*5*  j+®5 

is  a similar  4-point  average.  The  quantities  C n and  C n are 

yi1  Xi1 

approximated  as  follows: 

Cy"  -C]1  D+lD-i>]  + hJ<V">2 

3 I 1 


+ s.D  .R  ° + f.D  .Rn  + h (Rn  - R"  ) 
i ol  rz  i oj  rr  j rr  06  * 
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Cs”  -‘I1  [ -rj’n>o1  - <8im0l  - D«°-l  ] + »iDoiR° 


+ (f,D  . + h )Rn 

i oj  y rz  j 

The  finite-difference  approximation  to  the  divergence  of  the  mean 
velocity  field  at  n is  computed  from 

’’l't'i.W  " hrtfi-Bi(Dojuo 

+ WDoiH)U.3*  ' 

As  mentioned  in  Section  2.2  of  this  report,  we  use  the  NCAR  BLKTRI 
direct  Poisson  solver  to  solve  the  above  difference  equation  for 
Pn 

* 

D. 5 East  Boundary  Conditions 

During  the  first  half  time-step,  y derivatives  in  the  transformed 
and  "unsteady"  versions  of  Equations  3.10  and  3.11,  which  apply  at  the 
East  boundary,  are  differenced  implicitly,  and  x derivatives  in  the 
axial  convection  terms  are  approximated  with  first-order  backward  dif- 
ferences. The  following  equations,  which  apply  at  i * N , then  result: 
Vn+4  _ vn  ^ 

^t/2  ^ + ^Vol7^  + 8N1^JD-iVNj  + hjUNjVNj 

- Re_1[  + + wDoj  - hj]  c 

f jDoj  + 2hj  ^ * 


_ „n 

JU Jil  + f un  D + 8 w11  D w" 

At/2  jNjoj^Nj  SNWNj  -rVj 


m (pn  4,  p°  _ pn  pn  \ 

V N.j-Mi  + rN,J-»s  FN-J5,j-H5  N'Jj, j-Jj' 


+ Re'1  [flV-3  + <8J  + W*]  C - [£JD.3  + -3  ]R°.H1  • . 


I 


« 
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Bo  th  of  these  finite-difference  equations  can  be  expressed  in  the 
tridiagonal  form 


Vj^N.j-l  + V)Vj  + Sj^N.j+l  “ °Nj  * 

These  tridiagonal  systems  are  solved  for  and  W^5  after 

incorporation  of  the  boundary  conditions  at  j ■ 1 and  j » M . The 
transformed  version  of  Equation  2. 12  is  then  integrated  numerically  from  the 
axis  (since  U„  , ■ 0 ) to  obtain  . After  the  turbulence 

N,  1 Nj 

quantities  are  computed  at  time  level  n4*s  , Equation  3. 9 is  integrated 

_ iL 

numerically  for  PN  , which  is  used  as  the  East  boundary  condi- 

tion  in  the  calculation  of  . 

During  the  second  half  time-step,  x derivatives  are  differenced 
implicitly,  and  the  following  difference  equation  results  for 
with  a similar  one  for  w"*1  : 


n+i 

JLL 


_ Vn+Js 

ILL. 


n+l 


1/2 


+ f Un+JsD  V1*4*  + s Wtt*4sD  V 
j Nj  ojVNj  ®NWNj  -lVNj 


( 


c 


c 


c 


- Re 


-1 


(5v-j * s * -s  ] 
]' 


f .D  , + 2h,  R 
j oj  j re 


This  equation  can  be  rearranged  and  expressed  in  the  form 


.n+l  n+l  _ _ 

Vj^N-l.j  + BNj^Nj  °Nj 


for  each  value  of  j . The  application  of  this  kind  of  condition  to  the 
trldlagonal  equations  for  and  (i  < N)  is  straightforward. 

The  procedure  is  discussed  in  Roache  (1972).  Following  the  calculation 
of  and  (2  £ i £ N) , Equation  3-12  is  integrated  to  obtain 


Un+1 

NJ 

Un+1 

ij 


which  is  used  as  a boundary  condition  in  the  calculation  for 
(1  < N).  Finally,  after  the  calculation  of  the  turbulence 


C 
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quantities,  Equation  3-9  is  integrated  for 
boundary  condition  in  the  calculation  for 


pn+l 

N.j-Hs 

_n+l 

i+*S. 


which  is  used  as  a 
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